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ABSTRACT 


The  Telemetry  (TM)  Tracker  system  is  designed  to  determine  time,  space,  and 
position  information  (TSPI)  of  an  airborne  missile  by  detecting  its  telemetry  signal  at  a 
number  of  receiver  sites.  Doppler  frequency  measurements  are  converted  to  range 
differences  between  the  missile  and  receiver  sites,  whose  locations  are  known  in  three 
dimensions.  An  algorithm  then  utilizes  these  range  differences  to  obtain  missile  TSPI 
with  1 -meter  accuracy.  The  TM  Tracker  was  fielded  during  live  missile  firings  and 
measurements  indicated  that  the  desired  TSPI  accuracy  was  not  attained. 

This  thesis  examines  system  requirements  and  limitations  of  the  TM  Tracker  to 
obtain  TSPI  with  1-meter  accuracy.  The  theory  of  operation  and  components  of  the  TM 
Tracker  are  introduced.  Algorithms  used  in  computing  position  of  a  radiating  source 
from  range  differences  are  analyzed.  MATLAB  simulations  are  conducted  with  missile 
trajectory  data  to  determine  the  required  measurement  precision  and  signal-to-noise  ratio 
(SNR)  at  the  receiver  sites  to  obtain  1-meter  TSPI.  The  receivers’  45-degree,  3-dB  beam 
widths  are  then  implemented  to  observe  their  effects  on  TSPI  accuracy.  Simulations 
reveal  that  the  TM  Tracker  system  is  capable  of  producing  TSPI  with  1 -meter  accuracy 
provided  that  precise  frequency  measurements  and  adequate  SNR  values  are  available  at 
the  receiver  sites. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


vi 


TABLE  OF  CONTENTS 


I.  INTRODUCTION . 1 

H.  TELEMETRY  TRACKER  SYSTEM . 5 

A.  THEORY  OF  OPERATION . 5 

B.  SYSTEM  DESCRIPTION . 7 

C.  FUNCTIONAL  COMPONENTS . 9 

1 .  Control  and  Processing  Site . 9 

2.  Antennas . 

3.  IF  Receiver . H 

4.  GPS  Time  and  Frequency  Reference . 13 

5.  Postprocessing . 14 

D.  TM  TRACKER  FIELD  TESTING . 14 

IE.  SOURCE  ESTIMATION  USING  RANGE  DIFFERENCES . 19 

A.  INTRODUCTION . 19 

B.  FANG’S  ALGORITHM . 19 

C.  SMITH  -  ABEL  ALGORITHM . 23 

IV.  MATLAB  SIMULATIONS . 29 

A.  BACKGROUND . 29 

B.  PRELIMINARY  SIMULATIONS . 29 

C.  TM  TRACKER  SIMULATIONS . 34 

1 .  Missile  Trajectory  and  Remote  Receiver  Locations . 34 

2.  Computer  Algorithm  Implementation . 35 

a.  Trajectory  Data . 35 

b.  Missile  Velocity  Components . 35 

c.  TM  Tracker  Simulation  Block  Diagram . 37 

d.  Ranges  and  Initial  Range  Differences . 38 

e.  Doppler  Frequency  /  Differential  Range  Difference  Calculation . 38 

f.  Smith-Abel  Algorithm  Position  Estimation . 40 

g.  Range  Difference  and  Position  Error  Calculation . 40 

3 .  Algorithm  Results  Using  All  Receivers . 40 

4.  Algorithm  Results  Using  Available  Receivers . 45 

a.  Antenna  Beam  Width  Considerations . 45 

b.  Receiver  Time  Utilization . 45 

c.  Initial  Range  Differences  for  Added  Receivers . 47 

d.  Simulation  Results . 48 

V.  CONCLUSION . 55 

A.  SUMMARY . 55 

B.  FUTURE  RESEARCH . 56 

APPENDIX  A.  MATLAB  SOURCE  CODES . 57 

APPENDIX  B.  RECEIVER  AZIMUTH  AND  ELEVATION  ANGLES . 99 

LIST  OF  REFERENCES . 105 

INITIAL  DISTRIBUTION  LIST . 107 


vii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


LIST  OF  FIGURES 


1 .  Figure  2.1.  Intersecting  Hyperbolas  Using  Four  Receivers . 6 

2.  Figure  2.2.  TM  Tracker  Operation  Block  Diagram  From  Ref.  [1] . 8 

3.  Figure  2.3.  Control  and  Processing  Site  From  Ref.  [1] . 10 

4.  Figure  2.4.  TM  Tracker  Antennas  From  Ref.  [  1  ] . 11 

5.  Figure  2.5.  Antenna  Pattern  for  TM  Receiving  Antenna  From  Ref.  [1] . 12 

6.  Figure  2.6.  Single  Channel  Receiver  Block  Diagram  From  Ref.  [1] . 13 

7.  Figure  2.7.  Receiver  Locations  From  Ref.  [1] . 15 

8.  Figure  2.8.  Detailed  Primary  Receiver  and  Launcher  Locations  From  Ref.  [1] . 16 

9.  Figure  2.9.  TSPI  Accuracy  From  Ref.  [1] . 18 

10.  Figure  3.1.  Receiver  Stations  and  Estimated  Position . 20 

11.  Figure  3.2.  Notation  and  Geometric  Relations  for  ith  Receiver . 24 

12.  Figure  4. 1 .  Simple  Missile  Trajectory  with  Four  Receiver  Sites . 30 

13.  Figure  4.2.  Position  Error  vs.  SNR . 31 

14.  Figure  4.3.  X-Position  Error  vs.  RMS  Measurement  Error  . ; . 32 

15.  Figure  4.4.  T-Position  Error  vs.  RMS  Measurement  Error . 32 

16.  Figure  4.5.  Z-Position  Error  vs.  RMS  Measurement  Error . 33 

17.  Figure  4.6.  Total  Position  Error  vs.  RMS  Measurement  Error . 33 

1 8.  Figure  4.7.  Missile  Trajectory  and  Remote  Receiver  Locations . 34 

19.  Figure  4.8.  Missile  Velocity  Components  with  Roll-Rate  Modulation . 36 

20.  Figure  4.9.  Missile  Velocity  Components  with  Modulation  Removed . 36 

21 .  Figure  4. 10.  Block  Diagram  of  TM  Tracker  Simulation . 37 

22.  Figure  4.11.  Range  Difference  Error  vs.  Time  -  Ideal  Conditions . 41 

23.  Figure  4.12.  Position  Error  vs.  Time  -  Ideal  Conditions . 41 

24.  Figure  4.13.  Range  Difference  Error  vs.  Time  -  35  dB  SNR . 42 

25.  Figure  4.14.  Position  Error  vs.  Time  -  35  dB  SNR . 43 

26.  Figure  4. 15.  Range  Difference  Error  vs.  SNR . 44 

27.  Figure  4.16.  Position  Error  vs.  SNR . 44 

28.  Figure  4.17.  Azimuth  and  Elevation  Angles  for  Receiver  14 . 46 

29.  Figure  4.18.  Range  Difference  Error  vs.  Time  -  Ideal  Conditions . 48 

30.  Figure  4.19.  Position  Error  vs.  Time  -  Ideal  Conditions . 49 

31.  Figure  4.20.  Range  Difference  Error  vs.  Time  -  35  dB  SNR . 50 

32.  Figure  4.21 .  Position  Error  vs.  Time  -  35  dB  SNR . 50 

33.  Figure  4.22.  Range  Difference  Error  vs.  SNR . 51 

34.  Figure  4.23.  Position  Error  vs.  SNR . 52 

35.  Figure  B.l.  Reference  Receiver  Azimuth  and  Elevation  Angles . 99 

36.  Figure  B.2.  Receiver  lb  Azimuth  and  Elevation  Angles . 100 

37.  Figure  B.3.  Receiver  lp  Azimuth  and  Elevation  Angles . 100 

38.  Figure  B.4.  Receiver  10  Azimuth  and  Elevation  Angles . 101 

39.  Figure  B.5.  Receiver  11  Azimuth  and  Elevation  Angles . 101 

40.  Figure  B.6.  Receiver  12  Azimuth  and  Elevation  Angles . 102 

41.  Figure  B.7.  Receiver  14  Azimuth  and  Elevation  Angles . 102 

42.  Figure  B.8.  Receiver  15  Azimuth  and  Elevation  Angles . 103 

43.  Figure  B.9.  Receiver  16  Azimuth  and  Elevation  Angles . 103 


IX 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


x 


LIST  OF  TABLES 


1.  Table  2.1. 

2.  Table  4.1. 

3.  Table  4.2. 

4.  Table  A.  1. 


Receiver  Site  Antennas  and  Launch  Locations  From  Ref.  [1] 

Time  Utilization  of  Receivers  During  Missile  Trajectory . 

Summary  of  TM  Tracker  Simulations . 

MATLAB  Files . 


16 

47 

53 

57 


xi 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xii 


ACKNOWLDEGEMENT 


I  wish  to  express  my  sincere  gratitude  to  those  who  contributed  to  the  successful 
completion  of  this  project. 

I  would  like  to  thank  Greg  Velicer  and  the  personnel  at  the  Naval  Air  Warfare 
Center  -  Weapons  Division,  China  Lake,  CA  for  the  use  of  their  TM  Tracker  final  report 
and  the  expert  guidance  and  support  they  provided  throughout  the  past  nine  months.  I 
would  also  like  to  thank  Professor  David  Jenn  for  proofreading  my  thesis.  Additionally, 
many  thanks  go  out  to  my  thesis  advisor,  Professor  Curtis  Scheleher,  who  provided  me 
with  the  knowledge  and  guidance  to  program  the  simulations  in  MATLAB. 

Finally,  I  would  like  to  thank  my  wife  Robin  for  her  support  and  understanding 
during  the  many  long  hours  that  I  spent  at  the  campus  computer  lab  in  completing  this 
project. 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xiv 


I.  INTRODUCTION 


The  accurate  determination  of  time,  space,  and  position  information  (TSPI)  of  an 
airborne  vehicle  such  as  a  missile  has  been  a  valuable  asset  on  Department  of  Defense 
(DOD)  test  ranges  for  the  evaluation  and  testing  of  new  weapon  systems.  Four 
commonly  used  methods  in  determining  TSPI  include  high-speed  photography,  use  of  the 
Global  Positioning  System  (GPS),  laser  tracking,  and  radar  tracking.  While  the  spectrum 
of  TSPI  accuracy  ranges  from  inches  to  several  feet  with  the  use  of  these  methods,  they 
are  unable  to  meet  the  variance  of  field  test  requirements  over  all  DOD  test  ranges. 

High-speed  photography  uses  multiple  cameras  along  the  missile’s  trajectory  with 
high-speed  film  imagery.  These  cameras  are  coupled  with  a  tracking  mount  that  is  used 
by  an  operator  to  visually  track  the  airborne  missile.  TSPI  is  calculated  by  correlating  the 
angular  positions  of  the  cameras  with  the  missile’s  location  in  each  frame  of  film.  This 
method  can  produce  TSPI  data  with  several  feet  of  accuracy,  but  is  limited  in  the 
performance  of  its  operators  against  high-speed  vehicles  such  as  supersonic  missiles,  the 
capability  to  be  used  over  water,  and  its  use  with  multiple  targets.  GPS  uses  a 
constellation  of  satellites  to  determine  position,  but  is  limited  by  satellite  availability  at 
the  location  of  the  test  range.  Laser  tracking  is  able  to  produce  TSPI  with  inches  of 
accuracy  but  is  extremely  expensive  to  use  and  is  not  readily  portable  for  use  on  all  test 
ranges.  Radar  tracking  is  a  common  method  used  on  DOD  test  ranges  today.  Depending 
on  the  platform,  TSPI  accuracy  ranges  between  several  and  hundreds  of  feet.  Numerous 
missile  ranges  over  water  utilize  the  U.S.  Navy’s  SPY-1  phased-array  radar  for  TSPI 
calculation,  but  are  limited  by  its  availability  within  the  fleet  of  Navy  ships.  [Ref.  1] 
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The  Telemetry  (TM)  Tracker  system  was  developed  to  counter  the  deficiencies 
within  these  TSPI  collection  methods.  Its  theory  of  operation  allows  for  providing  TSPI 
data  in  a  wide  range  of  field  test  conditions,  including  the  ability  to  be  used  over  water 
and  land  ranges,  and  in  all  weather  conditions  that  allow  such  a  test  to  be  conducted.  The 
TM  Tracker  offers  portability  for  use  on  all  DOD  test  ranges  and  has  the  ability  to  track 
multiple  test  vehicles.  The  objective  in  its  development  and  design  was  TSPI  accuracy 
within  1  meter. 

The  TM  Tracker  determines  TSPI  of  an  airborne  missile  by  detecting  its  telemetry 
signal  at  a  number  of  remote  receiver  sites  whose  positions  are  known  in  three 
dimensions.  The  Doppler  shifted  signals  are  received  and  converted  to  range  differences 
between  the  remote  receiver  sites.  In  its  post  processing  system,  an  algorithm  utilizes 
these  range  differences  to  determine  the  telemetry  source  location  and  hence,  missile 
TSPI.  The  TM  Tracker  was  field  tested  at  White  Sands  Missile  Range  during  two  testing 
programs  in  August  and  September  1999  and  yielded  unfavorable  TSPI  data. 

This  thesis  examines  the  system  requirements  and  limitations  of  the  TM  Tracker 
to  obtain  TSPI  with  1 -meter  accuracy.  MATLAB  was  extensively  used  in  the 
simulations  of  the  TM  Tracker  system  along  with  missile  data  obtained  from  the  laser 
tracker  on  the  range. 

There  are  five  chapters  and  two  appendices  within  this  thesis.  In  Chapter  II,  the 
theory  of  operation  and  components  of  the  TM  Tracker  are  introduced  and  discussed. 
Range  difference  algorithms  used  in  simulations  for  determining  source  location  are 
analyzed  in  Chapter  IH  The  heart  of  this  project  is  contained  in  Chapter  IV.  Results  of 
MATLAB  simulations  of  source  location  algorithms  with  frequency  measurement  errors 
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and  noise  are  presented.  Simulations  of  the  operation  of  the  TM  Tracker  with  missile 
trajectory  data  and  all  remote  receivers  are  discussed.  These  simulations  are  then 
modified  to  implement  the  receivers’  45-degree,  3-dB  antenna  beam  widths,  which  limit 
their  ability  to  detect  the  telemetry  signal  along  the  missile’s  trajectory  and  ultimately,  the 
accuracy  of  TSPI.  Chapter  V  then  summarizes  and  concludes  this  work  and  discusses 
ongoing  and  future  research.  MATLAB  source  codes  used  for  simulations  are  listed  in 
Appendix  A.  Plots  of  azimuth  and  elevation  angles  vs.  time  for  the  nine  receivers  on  the 
missile  range  are  contained  in  Appendix  B. 
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II.  TELEMETRY  TRACKER  SYSTEM 


A.  THEORY  OF  OPERATION 

The  Telemetry  (TM)  Tracker  determines  time,  space  and  position  information 
(TSPI)  of  an  airborne  missile  equipped  with  a  telemetry  transmitter  by  measuring  range 
differences  with  respect  to  time  between  the  missile  and  the  receiver  sites  utilized  during 
missile  flight  whose  three-dimensional  positions  are  already  known  prior  to  launch. 
Range  differences  are  measured  between  receivers  utilized  along  the  missile’s  trajectory 
and  a  reference  receiver  within  the  constellation.  For  N  receivers  being  utilized,  there 
exists  N  -  1  range  differences  that  are  used  to  compute  missile  TSPI.  Range  differences 
are  computed  using  Doppler  frequencies  that  are  obtained  at  the  receiver  sites  along  with 
the  telemetry  signal  wavelength  and  sample  time.  By  utilizing  the  Doppler  difference 
between  receivers  and  the  initial  receiver  location,  the  range  difference  of  arrival 
(RDOA)  is  continuously  determined.  This  RDOA  is  then  used  within  an  algorithm  to 
determine  the  missile’s  position. 

The  range  difference  determined  between  two  receivers  restricts  the  source 
location  to  a  hyperbola  with  the  two  receivers  as  foci.  By  using  four  receiver  sites,  three 
independent  hyperbolas  are  generated  from  three  range  differences  and  their  intersection 
defines  the  source  location.  Figure  2.1  illustrates  three  intersecting  hyperbolas  that 
represent  the  missile’s  estimated  position  using  four  receiver  sites  in  the  same  plane. 
Each  of  the  three  hyperbolas  has  two  branches.  If  the  receivers  are  in  different  planes, 
the  estimated  position  is  defined  by  intersecting  hyperboloids. 
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Figure  2.1.  Intersecting  Hyperbolas  Using  Four  Receivers 

The  four  hyperbolas  in  Figure  2.1  show  a  favorable  geometric  configuration  of 
the  receivers  for  position  estimation  using  range  differences.  Their  intersection  results  in 
minimal  source  location  ambiguity  as  the  hyperbolas  become  orthogonal  to  one  another. 
The  geometric  configuration  of  the  receivers,  thus,  can  affect  the  amount  of  ambiguity  in 
determining  missile  position,  particularly  when  the  hyperbolas  are  mutually  parallel, 
representing  a  worst-case  scenario.  Noise  present  within  the  receivers  and  frequency 
measurement  inaccuracy  can  also  introduce  errors  in  the  estimated  position. 

An  incremental  approach  is  used  to  compute  receiver  range  tracks  for  later 
conversion  to  range  differences.  Initial  range  measurements  are  obtained  by  using  the 
missile’s  initial  location  and  the  surveyed  locations  of  each  receiver  site,  both  of  which 
are  already  known.  A  GPS  receiver  is  also  contained  at  each  receiver  site  for  a  time  and 
frequency  reference.  During  missile  flight  the  range  to  the  receivers  varies,  causing  a 
Doppler  frequency  shift  in  the  received  TM  signal  at  each  receiver.  The  Doppler 
frequency  shift  is  realized  at  each  receiver  by  the  difference  between  the  received  TM 
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signal  from  the  missile  and  the  frequency  reference  derived  from  GPS  for  each  time 
period.  Using  the  known  wavelength  of  the  TM  transmitter,  this  Doppler  shift  is 
converted  to  a  change  in  range.  Range  tracks  for  each  receiver  site  are  obtained  by 
integrating  the  changes  in  range  over  time  from  the  initial  ranges  obtained  before  launch. 
These  individual  range  tracks  are  converted  to  range  differences  with  respect  to  a 
reference  receiver  from  which  TSPI  is  obtained  using  an  algorithm  that  requires  at  least 
four  range  differences  from  five  receivers.  [Ref.  1]  This  algorithm  is  discussed  and 
analyzed  in  Chapter  III. 

Data  processing  is  conducted  upon  completion  of  missile  flight.  Data  from  each 
receiver  site  is  transmitted  to  the  control  and  processing  site  over  an  RF  modem  data  link. 
The  modulation  is  removed  from  the  data  to  obtain  the  carrier  frequency  with  Doppler 
shift  for  conversion  to  range  tracks  and  then  range  differences  to  obtain  missile  TSPI. 
[Ref.  1]  Overall  operation  of  the  TM  Tracker  is  shown  in  Figure  2.2. 

B.  SYSTEM  DESCRIPTION 

Remote  receiving  sites  and  a  control  and  processing  site  are  used  for  signal 
reception  and  data  processing  within  the  TM  Tracker  system.  Each  site  contains  an  RF 
modem  which  is  slaved  to  the  master  modem  at  the  control  and  processing  site  for  system 
timing  and  inter-site  communication.  The  control  and  processing  site  controls  the  TM 
Tracker  system  and  contains  a  Pentium-II  class  workstation  with  software  for  remotely 
controlling  receiver  sites  and  downloading  data  from  these  sites  for  processing  and 
conversion  to  missile  TSPI.  [Ref.  1] 
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Figure  2.2.  TM  Tracker  Operation  Block  Diagram  From  Ref.  [1] 

Prior  to  the  test,  the  TM  signal  receiver  within  each  receiving  site  is  tuned  to  the 
TM  frequency  of  2.3  GHz  with  a  command  from  the  control  and  processing  computer. 
The  received  TM  signal  is  heterodyned  to  an  intermediate  frequency  (IF)  of  10  MHz,  and 
drives  the  IF  phase  detector  in  the  receiver  while  preserving  any  Doppler  shift  within  the 
signal.  The  frequency  shift  keyed  (FSK)  modulation  of  the  TM  signal  is  also  measured 
by  the  receiver  with  a  “1”  bit  transmitted  when  the  frequency  is  higher  than  the  carrier, 
and  a  “0”  bit  transmitted  when  the  frequency  is  lower  than  the  carrier.  The  ratio  of  “1” 
bits  and  “0”  bits  is  used  to  remove  modulation  shifts  in  the  IF,  leaving  only  the  carrier. 
[Ref.  1] 
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All  receivers  within  the  TM  Tracker  system  must  have  time  and  frequency 
synchronization  to  prevent  errors  in  determining  the  position  of  the  missile.  Each 
receiver  site  contains  a  GPS  receiver  to  provide  a  10-MHz  frequency  reference  and  a  one 
pulse-per-second  time  reference  to  the  GPS  satellite  atomic  clocks  to  obtain  the  precision 
necessary  for  measuring  the  missile’s  Doppler  frequency.  Each  receiver  is  also  phase- 
locked  to  the  frequency  reference  to  maintain  frequency  synchronization  within  the 
system.  [Ref.  1] 

C.  FUNCTIONAL  COMPONENTS 

The  TM  Tracker  system’s  primary  functional  components  for  signal  reception  and 
conversion  to  TSPI  include  the  control  and  processing  site,  antennas,  IF  receiver,  GPS 
time  and  frequency  reference,  and  post  processing. 

1.  Control  and  Processing  Site 

The  control  and  processing  site  contains  a  computer  with  control  and  processing 
software,  an  uninterruptible  power  supply  (UPS),  and  an  RF  modem  for  communication 
with  remote  receiver  sites.  The  computer  controls  the  primary  settings  and  functions  of 
all  remote  receivers,  such  as  operating  frequency,  system  status,  and  data  collection.  It 
also  receives  and  processes  data  from  these  receivers  for  TSPI  generation.  High-speed 
computing  power  is  not  a  requirement  for  sufficient  operation  of  the  system.  A  300-MHz 
Pentium-II  desktop  and  a  233-MHz  Pentium-MMX  laptop  were  both  used  with 
satisfactory  system  operation.  [Ref.  1] 

The  control  and  processing  site’s  distance  from  the  remote  receivers  is  limited  by 
the  operating  range  of  the  RF  modem,  approximately  20  miles.  Line-of-sight  visibility 
between  the  RF  modem  antenna  and  each  receiver  site  is  required  for  proper  operation  of 
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the  communication  link.  The  link  operates  in  the  902  -  928  MHz  frequency  band  and  has 
a  115.2  kilobaud  data  rate.  [Ref.  1]  Figure  2.3  shows  the  control  and  processing  site 
utilized  with  a  laptop  computer. 


Figure  2.3.  Control  and  Processing  Site  From  Ref.  [1] 


2.  Antennas 

Four  different  antennas  are  used  within  the  TM  Tracker  system.  A  6-dB  gain,  S- 
Band  helical  antenna  with  right-hand  circular  polarization  is  used  for  TM  signal 
reception.  At  remote  receiver  sites,  a  10-dB  gain,  Yagi  antenna  is  used  for  RF  modem 
signal  reception  and  transmission.  A  5-dB  gain  vertical  monopole  antenna  is  utilized  at 
the  control  and  processing  site  (see  Figure  2.3)  for  the  RF  modem  link  because  of  its 
omni-directional  radiation  pattern,  allowing  communication  with  all  remote  receiver 
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sites.  The  GPS  receivers  at  each  remote  site  have  an  integrated  patch  antenna  and  low- 
noise  amplifier  combination.  [Ref.  1]  Figure  2.4  shows  the  antennas  used  at  each  remote 
receiver  site.  Figure  2.5  shows  the  antenna  pattern  for  the  S-Band  helical  antenna  for  TM 
signal  reception  with  a  3-dB  beam  width  of  approximately  45°. 


Figure  2.4.  TM  Tracker  Antennas  From  Ref.  [1] 


3.  IF  Receiver 

A  zero-crossing  counter  approach  of  the  received  frequency  modulated  (FM)  TM 
signal  is  used  within  the  IF  receiver.  Because  of  low  transmitter  power,  possible  terrain 
blockage  of  the  signal,  and  multipath  effects,  high  receiver  sensitivity  is  required.  A 
down-converter  method  is  utilized  with  an  independently  tunable  second  frequency 
conversion  for  an  IF  receiver  passband  of  1.5  MHz.  Using  two  separate  second  IF 
channels  allows  for  the  capability  to  track  two  different  transmitters  with  a  required 
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frequency  separation  of  13  MHz.  [Ref.  1]  Figure  2.6  shows  the  single  channel  block 
diagram  of  the  receiver. 


6dB  Gain  Helical  Antenna 
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Angle  (deg) 

Figure  2.5.  Antenna  Pattern  for  TM  Receiving  Antenna  From  Ref  [1] 

A  fixed  2.5-GHz  local  oscillator  (LO)  converts  the  2.2  -  2.3  GHz  TM  band  to  a 
200  -  300  MHz  first  IF  band.  A  second  LO  signal  is  provided  by  a  210  -  310  MHz 
synthesizer,  which  is  tunable  in  20-kHz  increments  and  forms  the  second  IF  channel 
centered  at  10  MHz.  The  GPS  10-MHz  reference  signal  is  used  by  all  LOs  to  provide  the 
required  synchronization  between  each  receiver.  Detailed  specifications  of  receiver 
components  can  be  found  in  Ref.  [1]. 
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Figure  2.6.  Single  Channel  Receiver  Block  Diagram  From  Ref.  [1] 

4.  GPS  Time  and  Frequency  Reference 

Four  10-MHz  frequency  reference  signals  and  a  communications  connection  are 
provided  by  the  TrueTime  Model  XL-AK  GPS  receiver.  When  power  is  first  applied  to 
the  unit,  it  determines  its  current  location  and  then  stabilizes  its  10-MHz  internal 
oscillator  to  the  GPS  system  clock.  Whether  tracking  one  or  six  satellites,  the  oscillator 
remains  close  to  its  stability  specification  of  5  x  10~n,  but  experiences  a  large  phase  drift 
of  50  cycles  when  a  receiver  changes  a  satellite  in  the  tracked  constellation.  During 
testing,  drifts  of  one  cycle  between  receivers  were  experienced  between  30  and  45 
seconds,  translating  into  a  frequency  offset  of  approximately  7.5  Hz  between  receivers. 
[Ref.  1] 
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5.  Post  Processing 

Post  processing  consists  of  converting  the  IF  receiver  data  to  range  data  and  then 
converting  this  range  data  to  TSPI.  Each  receiver’s  data  file  is  analyzed  to  remove 
duplicate  and  unreliable  data.  Duplicate  data  is  determined  by  analyzing  each  record’s 
time  stamp.  Unreliable  data  is  determined  by  the  received  signal  strength  indication 
(RSSI)  voltage  for  each  IF  count  within  each  receiver.  Incremental  counts  with  RSSI 
voltage  levels  below  a  specified  threshold  are  removed,  along  with  the  previous  and 
subsequent  counts.  Linear  interpolation  is  then  used  to  fill  in  the  missing  sample  points. 
[Ref.  1] 

Ranges  from  the  receiver  site  to  the  missile  are  obtained  by  integrating  the 
incremental  ranges  from  the  initial  range  from  launch  point  to  each  receiver,  requiring 
knowledge  of  the  launch  time.  Range  data  from  all  receiver  sites  are  then  combined  to 
compute  missile  TSPI.  Because  of  the  TM  transmitter  drift  error  within  all  range 
computations,  a  range  difference  algorithm  is  utilized  for  obtaining  the  missile’s  position, 
causing  these  drift  errors  to  cancel  out.  The  range  difference  algorithm  for  location 
estimation  developed  by  Julius  O.  Smith  and  Jonathan  S.  Abel,  explained  in  Chapter  III, 
is  used  to  compute  missile  TSPI. 

D.  TM  TRACKER  FIELD  TESTING 

The  TM  Tracker  was  field  tested  at  White  Sands  Missile  Range  during  two  testing 
programs  in  August  and  September  1999.  In  the  first  test  a  TM  transmitter  was  mounted 
on  a  cable  car  rolling  on  the  aerial  cable  of  the  Aerial  Cable  Range  (ACR).  Because  of 
the  transmitter  location,  several  receiver  sites  were  unable  to  receive  the  TM  signal 
resulting  in  unreliable  data.  [Ref.  1] 
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The  second  test  occurred  during  September  and  October  1999.  Missiles  with  TM 
sections  were  fired  against  targets  attached  to  the  cable  car  on  the  ACR.  Two  launch 
points  were  utilized,  with  approximate  ranges  of  3.0  km  and  4.7  km  from  the  target.  Nine 
receiver  sites  were  used  during  this  test.  Three  receiver  sites,  Rla,  Rib,  and  Rip,  were 
placed  near  the  launch  site  and  shared  the  same  receiver  site  enclosure.  The  other  six 
receivers  were  placed  at  other  positions  on  the  test  site.  [Ref.  1]  The  two-dimensional 
test-site  layout  and  the  layout  for  the  primary  receiver  (Rl)  and  3 -km  launch  point  are 
shown  in  Figures  2.7  and  2.8,  respectively.  Specific  coordinates  of  each  receiver  are 
listed  in  Table  2.1. 


Figure  2.7.  Receiver  Locations  From  Ref.  [1] 
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Launch  Point  1 
(4.7  km) 


Launch  Point  2 
(3.0  km) 


R1  a 


R1  b 


Rip 


RIO 


Rll 


R12 


R14 


R15 


R16 


X  (downrange)  feet 


-14863.3 


-9279.7 


-9243.4 


-9193.6 


-9278.7 


-8992.9 


1757.0 


-6458.5 


-3704.2 


607.0 


-5305.4 


Y  (crossrange)  feet 


-3676.5 


-3085.7 


-3085.7 


-3072.3 


-3201.8 


-2879.8 


430.4 


-2112.3 


-1214.7 


-3517.8 


2653.0 


Z  (elevation)  feel 


-250.2 


-122.3 


-127.0 


-125.0 


-92.1 


-109.3 


40.4 


-39.3 


-33.6 


208.3 


-153.6 


Table  2.1.  Receiver  Site  Antennas  and  Launch  Locations  From  Ref.  [1] 
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Receivers  near  the  launcher,  Rla,  Rib,  Rip,  and  RIO,  which  provide  the  initial 
missile  TSPI,  experienced  interference  with  the  received  TM  signal  around  launch  time. 
Several  sites  also  did  not  receive  a  strong  enough  signal  at  launch  time,  preventing  the 
determination  of  the  TM  receiver  frequency  offset  at  these  sites.  TSPI  was  only 
calculated  using  receiver  sites  Rip,  Rll,  R12,  and  R16,  after  the  interference  ended, 
while  the  laser  tracker  data  was  used  for  initial  missile  TSPI.  [Ref.  1] 

A  2-second  interval  of  data  with  minimal  receiver  dropouts  was  examined.  To 
compare  missile  TSPI  accuracy,  laser-tracking  data  of  the  missile  from  the  Sandia  Labs’ 
Laser  Tracker  was  utilized.  Figure  2.9  illustrates  the  total  position  error  for  x 
(downrange)  and  y  (crossrange)  dimensions  during  the  2-second  interval.  TSPI  data  of 
altitude  (z-dimension)  was  unreliable  and  is  not  shown.  [Ref.  1] 

The  results  of  the  live  fire  tests  were  far  from  optimal,  with  a  mean  position  error 
of  approximately  70  feet.  The  TM  Tracker’s  ability  to  obtain  1 -meter  position  accuracy 
requires  the  understanding  and  implementation  of  a  source  estimation  algorithm  within 
its  system  components.  Two  algorithms  that  utilize  range  differences  to  estimate  source 
location  are  discussed  and  analyzed  in  the  next  chapter,  one  of  which  is  used  within  the 
TM  Tracker’s  post  processing  system. 
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III.  SOURCE  ESTIMATION  USING  RANGE  DIFFERENCES 


A.  INTRODUCTION 

Position  estimation  of  a  radiating  source  by  analyzing  the  signals  received  at  a 
number  of  arbitrarily  spatially  separated  receiving  stations  using  the  range  differences  of 
the  propagating  signals  has  been  extensively  studied  [Refs.  2-5].  Signals  propagating 
from  the  source  arrive  at  the  receiving  stations  at  different  times  depending  on  the 
transmitter-to-receiver  geometry  and  the  medium  characteristics  through  which  the  signal 
is  propagating.  To  obtain  range,  signal  arrival  times  are  developed  at  all  receiver  stations 
relative  to  a  designated  reference  receiver  station.  These  time  difference  of  arrival 
(TDOA)  measurements  are  then  transformed  into  range  differences  between  receiver 
stations,  resulting  in  a  set  of  nonlinear  hyperbolic  equations  whose  solution  is  the 
estimated  position  of  the  radiating  source. 

Many  algorithms  have  been  developed  to  solve  the  nonlinear  hyperbolic  equations 
[Refs.  2-5].  Algorithms  developed  by  Bertrand  Fang  [Ref.  2]  and  Julius  Smith  and 
Jonathan  Abel  [Ref.  3]  were  used  in  MATLAB  simulations  (Chapter  IV)  and  are 
discussed  in  this  chapter.  Regardless  of  the  method  utilized,  reliable  TDOA 
measurements  mandate  precisely  synchronized  clocks  at  the  receiver  stations,  hence,  the 
use  of  GPS  receivers  in  the  TM  Tracker  system  [Ref.  1]. 

B.  FANG’S  ALGORITHM 

The  source  estimation  algorithm  developed  by  Fang  [Ref.  2]  was  utilized  during 
preliminary  feasibility  studies  for  1-meter  position  accuracy  of  the  TM  Tracker  system. 
A  system  with  three  receiver  stations  was  used,  with  the  origin  centered  on  one  of  the 
stations  as  shown  in  Figure  3.1. 
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From  Figure  3.1  the  range  differences  from  the  navigation  position  to  station  pairs  A,  B 


and  A,  C,  are  Rab  and  Rac,  respectively,  and  are  calculated  as 

V*2  +y2  +zl  +y2 +z2  =Kb 

(3.1) 

Jx2  +  y2+z2 -J(x-cJ  +(y-cyf  +  z2  =Rac 

(3.2) 

Squaring  and  simplifying  (3.1)  and  (3.2)  yields 

R Qb  “  b2  +  2 b  *  x  —  2 Rab x2  4-  y~  +  z2 

(3.3) 

Rl  -c2  +2cx*x  +  2cy*y  =  2Rac^jx2  +y2+z2 

(3.4) 

where  b  and  c  are  the  station  ranges,  c  =  «Jc2  +  c2  ,  and  *  denotes  multiplication. 
Equating  (3.3)  and  (3.4)  and  simplifying  results  in 


y=g*x+h  (3.5) 
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where 


8  = 


Kc*{blRab)-Cx 


2c„ 


(3.6) 


(3.7) 


Substituting  (3.5)  into  (3.3)  results  in 


where 


d*x2  +e*x  +  f  =  z2 

(3.8) 

d  =  4-(i>/Rj+g2} 

(3.9) 

e  =  f\-(b!  Rj}-2g*h 

(3.10) 

(3.11) 

The  navigation  position  vector  can  now  be  written  as  a  function  of  x  and  defines  a 
hyperbola  (d  >  0)  or  an  ellipse  (d  <  0)  with  mirror  symmetry  with  respect  to  the  receiver 
station  plane 

R  =  x*T +(g*x+h)j  ±-y/*/*;c2+e*x+ fk  (3.12) 


By  restricting  the  navigation  position  to  be  on  a  different  plane  than  the  base  stations 
results  in  the  position  becoming  the  solution  to  a  quadratic  equation.  For  position 
estimation  in  two  dimensions  x  and  y,  the  solution  to  (3.8)  for  x  is  found  by  applying  the 
Quadratic  Formula 


(3.12) 
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The  ambiguity  of  the  two  roots  in  (3.12)  can  be  resolved  if  knowledge  of  the  general 
location  is  known.  The  solution  fory  is  found  using  (3.5). 

Another  method  used  to  obtain  the  navigation  position  is  by  utilizing  a  fourth  base 
station,  Station  C'  (cx,  cy,  0),  which  results  in  a  navigation  position  that  is  at  the 
intersection  of  three  hyperboloids  or  ellipsoids  depending  on  the  value  of  d.  Computation 
follows  as  before  with  the  addition  of  another  range  difference  from  the  added  base 
station  to  the  reference  station  centered  at  the  origin 

-y/x^Z+z2  +z2  =R'ac  (3.13) 

Two  instances  of  (3.8)  result  for  Station  C  and  Station  C' 

d*x2  +e*x+f  =  z2  (3.14) 

d'*x2+e'*x  +  f'  =  z 2  (3.15) 

where  d',  e',  and /  are  computed  the  same  way  as  d,  e,  and /in  (3.9)  -  (3.11),  substituting 
R'ac,  cx',  and  cy'  for  Rac,  cx,  and  cy  for  calculating  g  and  h  in  (3.6)  and  (3.7).  Subtracting 
(3.15)  from  (3.14)  results  in  an  equation  in  terms  of  x 

(d -d')*x2  +(e-e')*x+(f-f')=  0  (3.16) 

Solving  for  x  results  in 

-(e-epl^-e')2  -4*  (d -/') 

2 (d-d')  K  '  ’ 

Solutions  fory  and  z  follow  from  (3.5)  and  (3.8). 

Simulations  using  Fang’s  algorithm  with  four  receiver  stations  are  discussed  in 

Chapter  IV.  While  Fang’s  algorithm  provides  an  exact  solution  to  the  hyperbolic 

equations,  it  cannot  take  advantage  of  utilizing  range  differences  from  additional 
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receivers.  A  maximum  of  three  range  differences  can  be  used  with  four  receivers.  The 
TM  Tracker,  however,  has  the  potential  of  eight  range  differences  available  from  nine 
receiver  sites.  Therefore,  a  range  difference  algorithm  that  can  use  all  available  range 
differences  in  determining  position  must  be  used. 

C.  SMITH -ABEL  ALGORITHM 

The  algorithm  for  source  location  developed  by  Julius  Smith  and  Jonathan  Abel 
[Ref.  3]  uses  a  spherical  interpolation  method  with  linear  least-squares  error 
minimization.  It  is  a  closed-form  method  that  uses  range  differences  referenced  to  a 
single  receiver  (reference  receiver)  with  N  -  1  range  differences  available  within  a 
constellation  of  N  receivers.  This  algorithm  is  used  within  the  TM  Tracker’s  post¬ 
processing  system. 

Let  N  denote  the  number  of  receivers  and  let  d#  denote  the  range  difference 
between  the  z'th  and  yth  receivers  and  the  source.  The  vectors  of  (x,  y,  z)  spatial 
coordinates  for  the  fth  receiver  and  the  source  are  denoted  xj  and  xj  respectively.  The 
distance  between  the  z'th  receiver  and  the  source  is  defined  by 

A=k--2*l  (318) 

and  the  distances  from  the  origin  to  the  points  x,  and  xj  are  i?,  and  Rs. 

The  reference  receiver  is  arbitrarily  chosen  as  Receiver  1  and  the  coordinate 
system  origin  is  translated  to  this  receiver  giving 

xj  =  0j=>  Rj  =  0  and  Dj  =  |x5 1  =  Rs  (3.19) 

Figure  3.2  illustrates  the  notations  and  geometric  relations  between  the  reference  receiver 
(Receiver  1),  the  z'th  receiver  (Receiver  z),  and  the  source.  The  distance  between  the  z'th 
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receiver  and  the  source,  Du  is  found  by  adding  the  range  difference  between  the  reference 
receiver  and  the  ith  receiver  to  distance  from  the  reference  receiver  to  the  source. 


RECEIVER  i 


SOURCE 


Figure  3.2.  Notation  and  Geometric  Relations  for  ith  Receiver 


Using  the  law  of  cosines,  A  can  be  expressed  as 

Df  =  (Rs  +  dn f  =  R?  +  R2 -2* |xf| * \xs | cosd  (3.20) 


where  0  is  the  angle  between  the  vectors  x*  and  xj.  The  angle  0  between  the  two  vectors 
Xj  and  Xj  is  defined  by 


cose  =  Iili-  (3.21) 

Substituting  (3.21)  into  (3.20)  and  cancelling  out  terms  results  in 

(Rs+dJ  =R-+R;  ~2xfxs  (3.22) 

Expanding  and  shifting  terms  to  the  right-hand  side  of  the  equation  yields 

s,  =  Rf  ~  d\  - 2 R,dn  - 2xf  x,  (3.23) 

where  s,  is  the  introduced  equation  error  due  to  imprecise  measurements  and  is 
minimized  in  a  least-square  sense  to  provide  an  estimate  of  the  true  source  position. 
Equation  (3.23)  can  be  written  in  matrix  form  for  N- 1  range  differences  as 
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where 
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Here  S  represents  the  matrix  of  receiver  spatial  coordinates  after  translation  of  the 
reference  receiver  to  the  origin.  The  least-squares  solution  for  &  given  Rs  is 

x,  =\s-wfe_-2R,d)  (3.25) 

where  *  denotes  complex  conjugation  and  complex  conjugate 

S'w  =  (sTs)~'sT  (3.26) 

yields  the  minimizer  of  sre.  The  range  differences  can  be  weighted  according  to  the 
relative  confidence  of  each  range  difference.  Weights  of  zero  can  be  used  when  certain 
range  differences  are  not  utilized.  This  technique  was  used  to  simulate  when  the  receiver 
does  not  have  the  missile  within  its  antenna  beam  (Chapter  IV).  This  weighted  equation 
error  is  erWe  and  is  minimized  for 

S’w  ={sTWsY  STW  (3.27) 

To  obtain  the  source  position  x*  by  minimizing  sTfVs,  Rs  must  be  allowed  to  vary  while 
maintaining  the  relation  Rs  =  |xs  |  because  the  range  from  the  reference  receiver  to  the 

source  is  not  typically  known  in  advance,  resulting  in  nonlinear  minimization.  This 
nonlinearity  can  be  approximately  eliminated  by  the  application  of  the  spherical 
interpolation  method. 


The  goal  of  the  spherical  interpolation  method  is  to  minimize  the  equation  error 
with  respect  to  Rs.  Substitution  of  (3.25)  into  (3.24)  results  in  a  linear  least-squares 
problem  for  finding  Rs  due  to  the  fact  that  the  least-squares  estimate  of  xs  given  Rs  in 
(3.25)  is  linear  in  Rs. 

e'  =  S-2Rsd-SS'rfe-2Rsd)  (3.28) 

where  s '  is  the  new  equation  error  which  is  linear  in  Rs.  Factoring  out  the  common  term 
6  -2Rsd  results  in 

e'  =  (l-SS^-2Rsd)  (3.29) 

where  /  is  the  identity  matrix. 

It  follows  from  (3.27)  that 

SS;  =  s(srWS)“'srW  =  Ps  (3.30) 

and 

P^I-P^I-SS;  (3.31) 

where  SS^  defines  the  TV  -  7  by  TV  -  7  symmetric  matrices,  and  Ps  andP/  are  rank  3 

projection  matrices  that  remove  both  orthogonal  components  and  components  within  the 
space  spanned  by  the  columns  of  S,  respectively. 

If  four  receivers  are  utilized,  three  range  differences  are  produced  and  Ps  =  /.  The 
equation  error  s' becomes  zero  for  all  values  of  Rs  resulting  in  ambiguity  for  the  source 
location.  Because  of  this  limitation,  a  minimum  of  five  receivers  producing  four  range 
differences  must  be  used  for  this  method  to  be  effective  in  non-linear  minimization. 

Assuming  the  use  of  at  least  five  receivers  and  four  range  differences,  substitution 
of  (3.31)  into  (3.29)  results  in 
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e'  =  Ps1(8_-2Rsd) 


(3.32) 


Utilizing  the  weighting  matrix  W  produces 

s  ,rWs '  =  (£  -  2 Rs  d)T  P^WP f  (5  -  2 Rs  d)  (3.33) 


A  form  of  weighted  least  squares  results  when  minimizing  (3.33)  with  respect  to  Rs 
where  the  weighting  matrix  P^WP^  is  of  rank  N  -  4.  The  three  missing  dimensions 

within  this  matrix  result  from  the  degrees  of  freedom  removed  by  substituting  (3.25)  in 
(3.24)  for  the  spatial  coordinates  of  the  source  Xg.  The  minimizing  value  of  Rs  for  (3.33) 
is 


~  _  dTP£WP& 
1  2drPs±WPsLd 


(3.34) 


Substituting  (3.34)  into  (3.25)  results  in  the  source  location  estimate  as 

x,  =±S-w{e-2R,d)=±(sTWSySTw(l-2R,d)  (3.35) 

The  source’s  coordinates  are  based  on  a  coordinate  system  with  the  reference  receiver  at 
the  origin.  To  obtain  the  source’s  true  spatial  coordinates,  the  coordinate  system  must  be 
shifted  back  by  translation. 

The  Smith-Abel  algorithm  for  source  location  estimation  was  utilized  in  all 
MATLAB  simulations  involving  the  TM  Tracker.  It  was  programmed  into  MATLAB 
and  used  as  the  function  smith  abel.m.  The  MATLAB  source  code  for  this  function  is 
listed  in  Appendix  A. 
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IV.  MATLAB  SIMULATIONS 


A.  BACKGROUND 

Several  questions  were  addressed  in  evaluating  the  potential  for  achieving  1-meter 
TSPI  accuracy  with  the  TM  Tracker  in  its  current  configuration,  including: 

1.  What  range  difference  measurement  accuracy  is  required  at  each  remote 
receiver  site  to  provide  1-meter  missile  TSPI  accuracy? 

2.  What  signal-to-noise  ratio  (SNR)  is  required  at  the  remote  receiver  sites  to 
obtain  1 -meter  missile  TSPI  accuracy  when  all  receivers  are  used  during  the 
entire  missile  trajectory? 

3.  What  are  the  effects  of  the  receivers’  antenna  beam  widths  on  TSPI  accuracy? 

4.  How  can  the  receivers  be  utilized  so  that  during  the  entire  missile  trajectory  at 
least  five  receivers  are  able  to  receive  the  telemetry  signal? 

5.  What  SNR  is  required  at  the  receivers  to  achieve  1 -meter  TSPI  when  only 
those  that  can  receive  the  telemetry  signal  within  their  respective  3-dB  beam 
widths  are  utilized? 

These  questions  were  answered  by  modeling  and  conducting  simulations  of  the  TM 
Tracker  system.  MATLAB  source  codes  for  all  simulations  are  listed  in  Appendix  A. 

B.  PRELIMINARY  SIMULATIONS 

Simulations,  in  MATLAB  using  Fang’s  algorithm  [Ref.  2]  were  conducted  to 
evaluate  the  required  measurement  accuracy  and  limitations  of  the  TM  Tracker  system. 
The  initial  approach  utilized  four  receiver  sites  in  the  X-Y  plane  that  were  placed  behind 
the  missile  trajectory  for  continuous  TM  signal  reception  throughout  missile  flight.  A 
simple  missile  trajectory  consisting  of  a  Mach  1  missile  flying  at  a  constant  altitude  of  1 
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km  was  used.  Figure  4.1  shows  the  receiver  locations  and  the  missile  trajectory  used  in 
preliminary  simulations. 


z 


Figure  4.1.  Simple  Missile  Trajectory  with  Four  Receiver  Sites 

TDOA  and  resulting  range  differences  were  obtained  using  the  frequency 
difference  of  arrival  (FDOA)  method  by  Chestnut.  [Ref.  5]  The  theoretical  lower  bound 
on  Doppler  frequency  measurement  was  used,  described  by  the  following  equation: 


n 


V  SNR 


7t-T 


(4.1) 


where  SNR  is  the  signal-to-noise  ratio  of  the  receiver  at  which  the  frequency 
measurement  is  made  and  T  is  the  sampling  period.  [Ref.  6]  The  theoretical  position 
accuracy  (mean  error)  for  a  given  value  of  SNR  was  computed  using  this  model.  Monte 
Carlo  simulation  was  utilized  due  to  the  random  noise  added  to  the  TDOA  measurements 
within  the  algorithm,  simulating  real  world  conditions  (fdoafang4nmc.m ).  Figure  4.2 
shows  the  mean  error  for  SNR  values  from  15  -  40  dB  along  with  a  curve  fit  of  the  data. 
The  results  from  this  model  established  the  feasibility  of  obtaining  1 -meter  position 
accuracy  using  TDOA  and  range  differences,  provided  that  sufficient  SNR  values  are 
obtained  at  each  receiver. 
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Position  Error  vs.  SNR 


SNR  -  dB 

Figure  4.2.  Position  Error  vs.  SNR 


The  effects  of  TDOA  measurement  accuracy  on  position  accuracy  were  also 
computed  using  Fang’s  algorithm  and  the  receiver  site  configuration  in  Figure  4.1.  A 
Monte  Carlo  simulation  using  10,000  iterations  for  each  value  of  root-mean-squared 
(RMS)  TDOA  measurement  error  was  used  to  compute  ^-position,  y-position,  z-position, 
and  total  position  error  ( tdoafang4nmc.m ),  illustrated  in  Figures  4.3  -  4.6,  respectively. 
Introduced  errors  in  TDOA  produced  range  difference  inaccuracies  which,  when  utilized 
in  Fang’s  algorithm,  resulted  in  increased  errors  in  position. 
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As  shown  in  the  previous  figures,  TDOA  measurement  accuracy  must  be  on  the 
order  of  1.2  inches  to  provide  1 -meter  position  accuracy,  with  the  z-position  accuracy 
being  the  most  critical  measurement  for  this  trajectory.  The  Fang  algorithm,  however, 
can  only  be  used  with  a  four-receiver  configuration  and  does  not  allow  the  addition  of 
more  receivers  for  possible  improved  TSPI  accuracy.  Because  of  this  limitation,  the 
Smith  -  Abel  algorithm  was  utilized  for  all  subsequent  simulations  of  the  TM  Tracker, 
allowing  for  the  use  of  nine  available  receiver  sites  that  were  on  the  missile  range. 

C.  TM  TRACKER  SIMULATIONS 

1.  Missile  Trajectory  and  Remote  Receiver  Locations 
The  remote  receiver  configuration  and  the  missile  trajectory  used  during  TM 
Tracker  testing  at  White  Sands  Missile  Range  were  utilized  for  all  simulations  and  are 
depicted  in  Figure  4.7.  Receivers  are  numbered  corresponding  to  Figure  2.8. 


Missile  Trajectory  and  Remote  Receiver  Locations 


Figure  4.7.  Missile  Trajectory  and  Remote  Receiver  Locations 
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2.  Computer  Algorithm  Implementation 

a.  Trajectory  Data 

The  heart  of  the  project  was  to  simulate  operation  of  the  TM  Tracker  using 
missile  trajectory  data  provided  from  the  laser  tracker  on  the  missile  range.  Spatial 
coordinates  and  time  data  were  resolved  by  the  laser  tracker  at  intervals  of  0.01  sec.  The 
trajectory  data  from  the  laser  tracker  was  stored  in  a  4  x  681  matrix.  The  four  rows 
within  this  matrix  represent  the  time  t  in  seconds  and  the  missile’s  spatial  coordinates  x, 
y,  and  z,  in  feet.  The  681  columns  represent  time  and  position  data  of  the  missile,  with 
the  first  and  last  columns  representing  t  =  0.01  sec  and  t  =  6.81  sec,  respectively.  Data  up 
through  t  =  6.42  sec  was  used  in  all  TM  Tracker  simulations  due  to  ambiguities 
encountered  in  the  data  near  the  end  of  the  missile’s  trajectory. 

b.  Missile  Velocity  Components 

The  trajectory  data  provided  by  the  laser  tracker  on  the  White  Sands 
Missile  Range  was  processed  through  a  SIMULINK  model  and  revealed  a  16-Hz 
modulation  component  in  the  x,  y,  and  z-component  velocities,  relating  to  the  missile’s 
roll  rate.  Because  Doppler  differential  measurements  were  to  be  used  for  obtaining 
TDOA  and  range  differences,  the  velocity  components  had  to  be  smoothed  to  eliminate 
missile  roll-rate  modulation  effects.  Figures  4.8  and  4.9  show  missile  velocity 
components  before  and  after  a  14th  degree  polynomial  least-mean-squares  fit  was  applied 
to  the  trajectory  data.  Velocity  components  were  found  by  computing  the  derivatives  of 
the  three  equations  for  missile  position  with  respect  to  time 
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Missile  Velocity  Components 


Figure  4.8.  Missile  Velocity  Components  with  Roll-Rate  Modulation 


Missile  Velocity  Components 


Figure  4.9.  Missile  Velocity  Components  with  Modulation  Removed 


c.  TM  Tracker  Simulation  Block  Diagram 

The  flow  of  calculations  for  simulation  of  the  TM  Tracker  in  MATLAB  is 
shown  in  the  block  diagram  of  Figure  4.10.  Subsequent  sections  discuss  the  calculations 
depicted  in  this  figure. 


Figure  4.10.  Block  Diagram  of  TM  Tracker  Simulation 
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d.  Ranges  and  Initial  Range  Differences 

The  ranges  from  each  receiver  to  the  reference  are  computed  before 
missile  flight  begins  because  these  values  remain  constant  and  are  one  of  the  inputs  to  the 
Smith-Abel  algorithm.  Range  differences  are  found  using  an  incremental  integrating 
approach.  The  initial  range  differences  are  calculated  prior  to  missile  launch  to  establish 
a  base  to  which  average  differential  range  differences  are  added  at  each  time  interval. 

e.  Doppler  Frequency  / Differential  Range  Difference  Calculation 

A  phase-detector  approach  for  frequency  measurement  was  used  for  all 

TM  Tracker  simulations  along  with  the  theoretical  lower  bound  on  frequency 
measurement  error  due  to  noise  (see  Equation  4.1)  because  of  the  inherent  limitations  of 
the  zero-crossing  counter  method  for  determining  frequency.  This  method,  implemented 
within  the  TM  Tracker,  requires  that  the  receivers  must  be  in  full-view  of  the  missile 
throughout  the  entire  trajectory  because  the  zero-crossing  count  begins  at  time  zero.  This 
also  introduces  a  uniform  phase  error  over  ±  360  degrees,  equivalent  to  0.130  m  error  in 
range.  Even  with  no  phase  measurement  noise  within  the  system,  position  error  due  to 
the  zero-crossing  counter  method  alone  was  on  the  order  of  10  feet.  This  error  only 
became  worse  when  phase  measurement  noise  and  receiver  noise  were  added  to  the 
system.  [Ref.  7] 

The  Doppler  frequency  vector  fd  detected  by  a  stationary  receiver  from  a 
source  with  a  wavelength  A  and  a  velocity  vector  v  is  defined  as 

H  (4-2) 
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To  compute  the  Doppler  frequencies  at  each  receiver,  fd  nc ,  the  scalar  component  of  fd 
in  the  direction  of  each  receiver  is  found  by  using  the  dot  product 


f  =/•  — 
fdrec  1  \§\ 


(4.3) 


where  B  is  the  vector  from  the  source  to  the  receiver  and  \B\  is  the  vector  magnitude. 


Each  scalar  Doppler  frequency  is  converted  to  a  differential  range  difference  (RD).  The 
RD  between  the  reference  and  ith  receivers  is  defined  by 

RD,=-A-T-{fdl-fdi)  (4.4) 

where  A  is  the  source  wavelength,  T  is  the  sample  time  (0.01  sec),  and^;  and  fa  are  the 
respective  Doppler  frequencies  at  the  reference  receiver  and  ith  receiver.  The  negative 
sign  ensures  a  positive  RD  value;  (/rfl  -fdi)  will  always  be  negative  because  the  missile 
is  moving  away  from  the  reference  receiver  through  its  entire  trajectory,  producing  lower 
relative  Doppler  frequencies  at  the  reference  receiver. 

Two  differential  RDs  are  then  computed  for  each  receiver  used.  The  first 
RD  is  based  on  missile  position  and  velocity  at  the  present  sample  time.  The  second  RD 
is  based  on  missile  position  and  velocity  at  the  next  sample  time  (0.01  sec  later).  These 
sets  of  RDs  are  used  to  determine  average  differential  RDs  which  are  added  to  the  base 
RDs  determined  previously.  The  base  RDs  are  then  reset  to  these  values  and  the  process 
is  repeated  through  the  missile’s  trajectory.  Geometric  RDs  are  also  calculated  to 
determine  the  error  in  RD  computation  using  missile  position  data  for  range  calculations. 
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f  Smith- Abel  Algorithm  Position  Estimation 

The  missile’s  estimated  position  at  the  next  time  interval  is  determined 


using  the  Smith- Abel  algorithm  ( smithjabel.m, ).  Inputs  to  this  algorithm  include  receiver 
location  coordinates,  RDs,  ranges  from  receivers  to  the  reference  receiver,  and  weight 
values.  For  all  simulations,  weights  were  only  assigned  values  of  one  or  zero  for  receiver 
RDs  that  were  used  or  not  used  during  the  missile’s  trajectory. 

g.  Range  Difference  and  Position  Error  Calculation 
The  range  difference  eixor  is  determined  by  computing  the  magnitude  of 
the  vector  from  the  geometric  RD  to  the  RD  obtained  using  incremental  RDs  from 
Doppler  frequencies.  Position  error  is  defined  as  the  magnitude  of  the  vector  from  the 
Smith-Abel  position  estimation  to  the  missile  position  coordinates  at  the  next  time 
interval.  These  errors  are  computed  for  each  sample  point  along  the  missile’s  trajectory 
and  stored.  Mean  values  for  these  errors  are  then  calculated  at  the  end  of  missile  flight. 

3.  Algorithm  Results  Using  All  Receivers 

Simulations  of  the  TM  Tracker  were  first  conducted  using  all  nine  receivers 
available  on  the  missile  range  ( trajall.m ).  Ideal  conditions  with  the  absence  of  noise  were 
implemented  in  these  simulations.  Because  there  was  no  noise  present  within  the  system, 
any  value  of  SNR  used  within  the  program  yielded  the  same  results.  Figures  4.11  and 
4.12  show  the  range  difference  and  position  errors  that  occurred  during  the  missile’s 
trajectory. 
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Range  Difference  Errors.  Time 


3  4 

Time  -  sec 


Figure  4.11.  Range  Difference  Error  vs.  Time  -  Ideal  Conditions 

Position  Error  vs.  Time 


Time  -  sec 

Figure  4.12.  Position  Error  vs.  Time  -  Ideal  Conditions 


Mean  range  difference  and  position  errors  were  quite  low,  0.046  ft  and  0.168  ft, 
respectively.  These  errors  represent  the  best  obtainable  accuracy  of  the  TM  Tracker 
under  perfect  conditions.  Real-world  conditions  with  random  white  Gaussian  noise 
present  within  the  Doppler  frequency  measurements  were  simulated  to  determine  effects 
on  range  difference  and  position  errors  ( trajalln.m ).  The  resulting  errors  with  receiver 
SNR  of  35  dB  are  shown  in  Figures  4.13  and  4.14. 


Range  Difference  Error  vs.  Time 
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Position  Error  vs.  Time 


Figure  4.14.  Position  Error  vs.  Time  -  35  dB  SNR 


The  values  were  a  0.123-ft  mean  error  in  range  difference  and  a  1.975-ft  mean 
error  in  position.  Because  of  the  randomness  of  the  noise,  these  errors  would  vary  with 
each  subsequent  run  of  the  program  using  the  same  SNR  value.  Monte  Carlo  simulation 
with  600  iterations  for  each  value  of  SNR  from  15  -  40  dB  was  conducted  to  determine 
the  range  difference  and  position  errors  over  a  wider  sample  of  data  ( trajallnmc.m ).  The 
results  are  shown  in  Figures  4.15  and  4.16. 
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Range  Difference  Error  vs.  SNR 
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Figure  4.15.  Range  Difference  Error  vs.  SNR 


The  results  from  these  figures  indicate  that  the  remote  receivers  require  an  SNR 
value  of  approximately  39  -  40  dB  to  obtain  TSPI  accuracy  within  1  meter  (about  3.28 
ft).  These  simulations,  however,  assume  that  all  receivers  can  detect  the  missile’s 
telemetry  signal  throughout  its  entire  trajectory.  Because  of  the  receivers’  locations  on 
the  range  and  their  antenna  beam  pattern  limitations,  this  is  not  possible. 

4.  Algorithm  Results  Using  Available  Receivers 

a.  Antenna  Beam  Width  Considerations 

All  receivers  have  a  3-dB  beam  width  of  approximately  45°  (see  Figure 
2.6).  To  determine  each  receiver’s  ability  to  obtain  the  telemetry  signal,  the  missile’s 
trajectory  had  to  be  converted  to  spherical  coordinates  and  referenced  to  each  receiver  in 
terms  of  azimuth  and  elevation.  Range  was  not  used  because  the  trajectory  did  not  pose 
any  range  limitations  on  the  receivers.  MATLAB  was  used  to  calculate  values  for 
azimuth  (theta)  and  elevation  (phi)  for  each  receiver  during  the  missile’s  trajectory 
(i angles.m ).  Because  MATLAB’s  conversion  from  Cartesian  to  spherical  coordinates  is 
referenced  to  the  origin,  the  coordinate  system  was  translated  such  that  the  origin  was 
centered  on  each  receiver,  resulting  in  also  translating  the  missile’s  trajectory  so  that  the 
relative  geometry  between  the  receiver  and  the  missile  was  maintained. 

b.  Receiver  Time  Utilization 

Plots  of  azimuth  and  elevation  versus  time  were  then  used  to  determine 
what  period  of  time  the  receivers  could  receive  the  telemetry  signal  in  their  respective  3- 
dB  beam  widths.  Since  the  antennas  are  stationary,  only  one  time  period  could  be 
utilized  for  each  receiver.  This  time  period  is  limited  by  a  45°  change  in  either  azimuth 
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or  elevation,  whichever  occurs  first.  Figure  4.17  shows  a  plot  of  these  angles  vs. 
trajectory  time  for  Receiver  14  and  how  its  time  period  was  determined. 


Angles  for  Receiver  14 


Time  -  sec 


Figure  4.17,  Azimuth  and  Elevation  Angles  for  Receiver  14 

In  the  case  of  Receiver  14,  the  time  period  was  chosen  at  the  start  of  the  missile’s 


trajectory  and  the  elevation  angle  phi  was  the  first  to  reach  its  45°  deviation  limit  at  about 
3.29  sec.  Four  receivers,  la,  lb,  lp,  and  10,  were  able  to  retain  the  telemetry  signal 
through  the  entire  trajectory  because  azimuth  and  elevation  angles  never  changed  by 
more  than  45°  during  missile  flight.  Many  different  scenarios  were  conducted  for  the 
other  five  receivers  to  determine  the  best  way  to  implement  their  respective  time  periods 
within  the  algorithm  to  achieve  the  lowest  range  difference  and  position  errors.  Plots  of 
azimuth  and  elevation  angles  vs.  time  for  each  receiver  are  listed  in  Appendix  B.  The 
time  period  utilization  of  receivers  for  subsequent  simulations  is  shown  in  Table  4.1. 
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Time 

Receivers  Used  or  Lost 

Start -2. 11  sec 

Use  Receivers  la,  lb,  lp,  10, 11, 12,  14, 15 

2.12  sec 

Lose  Receiver  12 

2.12  sec  -  3.28  sec 

Use  Receivers  la,  lb,  lp,  10,  11,  14,  15 

3.29  sec 

Lose  Receiver  14 

3.29  sec  -  4.48  sec 

Use  Receivers  la,  lb,  lp,  10, 11, 15 

4.49  sec 

Add  Receiver  16,  Lose  Receiver  15 

4.49  sec  -5.94  sec 

Use  Receivers  la,  lb,  lp,  10, 11, 16 

5.95  sec 

Lose  Receiver  1 1 

5.95  sec  -  End 

Use  Receivers  la,  lb,  lp,  10,  16 

Table  4.1.  Time  Utilization  of  Receivers  During  Missile  Trajectory 
c.  Initial  Range  Differences  for  Added  Receivers 

The  weights  within  the  Smith- Abel  algorithm  allow  for  changing  the  range 
differences  used  to  estimate  missile  position  during  its  trajectory.  Values  of  one  are 
applied  to  range  differences  in  use,  values  of  zero  are  applied  to  range  differences  not 
used  for  position  estimation.  Removing  receivers  during  the  trajectory  poses  no 
problems  in  the  simulation  as  long  as  a  minimum  of  five  receivers  (four  range 
differences)  is  maintained.  Adding  receivers  after  time  start,  however,  requires  special 
consideration  in  terms  of  computing  the  initial  base  range  differences  to  which  the 
incremental  range  differences  are  added  in  subsequent  time  intervals. 

When  a  receiver  is  added  after  missile  flight  has  begun,  an  initial  base 
range  difference  must  be  calculated  at  that  instant  based  on  the  missile’s  estimated 
position  determined  from  the  Smith-Abel  algorithm.  This  is  due  to  the  fact  that  an  initial 
range  difference  was  not  computed  prior  to  missile  launch  for  this  receiver  because  it 
could  not  receive  the  telemetry  signal  at  that  time.  Incremental  range  differences  are  then 
added  to  this  base  as  in  the  previous  simulations.  Because  this  base  is  derived  from  an 
estimated  position  that  contains  some  error,  position  errors  after  this  added  receiver  will 
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Range  Difference  Error  -  ft 


become  progressively  worse,  therefore,  it  is  best  to  add  receivers  as  late  in  the  trajectory 
as  feasible  to  prevent  the  rapid  accumulation  of  error  during  missile  flight. 


d.  Simulation  Results 

Receivers  were  utilized  in  MATLAB  as  shown  in  Table  4.1  in  ideal 
conditions  with  the  absence  of  noise  ( trajgroup.m ).  Figures  4.18  and  4.19  show  the  range 
difference  and  position  errors  that  occurred  during  missile  flight. 


Range  Difference  Error  vs.  Time 


Figure  4.18.  Range  Difference  Error  vs.  Time  -  Ideal  Conditions 
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Position  Error  vs.  Time 


Figure  4.19.  Position  Error  vs.  Time  -  Ideal  Conditions 

Mean  values  for  range  difference  and  position  error  were  0.054  ft  and 
2.200  ft,  respectively.  A  jump  in  range  difference  error  can  be  noticed  in  Figure  4.18  at 
4.49  sec  when  Receiver  16  was  added  due  to  the  base  range  difference  being  calculated 
from  the  missile’s  estimated  position  at  that  point.  Range  difference  error  then  remains 
stable  from  that  point  until  the  end  of  the  trajectory.  Figure  4.19  shows  slight  increases  in 
positional  error  as  receivers  were  lost  during  the  trajectory,  but  a  drastic  increase 
occurred  when  Receiver  11  dropped  out  at  5.95  sec.  Large  increases  in  position  error 
also  occurred  at  the  end  of  trajectory  during  testing  of  other  receiver  utilization  methods. 
White  Gaussian  noise  was  added  to  Doppler  frequency  measurements  in  this  simulation 
to  represent  real-world  conditions  {trajgroupn.m).  Figures  4.20  and  4.21  show  resulting 
range  difference  and  position  errors  with  receiver  SNR  of  35  dB. 
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Range  Difference  Error  -  ft 


Range  Difference  Error  vs.  Time 
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Figure  4.20.  Range  Difference  Error  vs.  Time  -  35  dB  SNR 

Position  Error  vs.  Time 


Time  -  sec 

Figure  4.21.  Position  Error  vs.  Time  -  35  dB  SNR 


Added  noise  caused  a  slight  increase  in  mean  range  difference  error,  0.169 
ft,  and  a  larger  increase  in  mean  position  error,  7.681  ft.  Even  with  the  randomness  of  the 
noise,  large  increases  in  position  error  tend  to  occur  after  approximately  4.5  sec.  Monte 
Carlo  simulation  was  conducted  with  600  iterations  for  each  value  of  SNR  from  15-40 
dB  to  better  estimate  the  effects  using  a  larger  sample  of  data  (trajgroupnmc.m).  Figures 
4.22  and  4.23  depict  the  errors  that  occurred. 


Range  Difference  Error  vs.  SNR 


Figure  4.22.  Range  Difference  Error  vs.  SNR 
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Position  Error  vs.  SNR 


Figure  4.23.  Position  Error  vs.  SNR 

Figures  4.22  and  4.23  show  that  receiver  SNR  values  need  to  be  above  40 
dB  to  obtain  1-meter  TSPI  accuracy  for  the  method  used  in  Table  4.1.  It  is  important  to 
note  that  an  enormous  amount  of  methods  to  utilize  the  receivers  during  the  missile’s 
trajectory  exists.  The  receiver  time  windows  used  in  these  simulations  were  chosen  after 
many  different  combinations  of  time  windows  were  tested.  During  the  testing  of  these 
combinations,  position  errors  increased  when  receivers  were  added  near  the  beginning  of 
the  missile’s  trajectory  due  to  the  rapid  accumulation  of  error  that  occurs  from  estimating 
the  initial  range  difference  of  the  added  receiver.  If  this  accumulation  is  allowed  to  start 
early  in  the  missile’s  trajectory,  a  larger  error  will  result  by  the  time  the  end  of  the 
trajectory  is  reached.  Adding  receivers  later  in  the  trajectory  minimizes  this  error,  but  at 
the  same  time,  a  minimum  of  five  receivers  must  always  be  maintained  which  could  limit 
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the  use  of  this  method.  The  best  implementation  with  the  least  amount  of  error  is  when 
all  nine  receivers  are  utilized  throughout  the  entire  trajectory  as  shown  in  Figures  4.15 
and  4.16.  Errors  in  computing  added  receivers’  initial  range  differences  are  eliminated 
because  all  receivers  are  used  from  time  start  until  the  end  of  the  trajectory.  Relocating 
receivers  along  with  the  geographic  conditions  on  the  range,  however,  may  preclude  the 
use  of  this  option. 

A  table  was  developed  to  summarize  the  results  of  TM  Tracker 
simulations.  Table  4.2  lists  the  simulation  description  and  resulting  mean  range 
differences  and  position  errors  that  were  produced. 


TM  Tracker 
Simulation  Description 

MATLAB  File 
Name 

Mean  Range 
Difference  Error  (ft) 

Mean  Position 
Error  (ft) 

All  receivers  used  with  no 
noise 

trajall.m 

0.046 

0.168 

All  receivers  used  with 
noise,  35  dB  SNR 

trajalln.m 

0.123 

1.975 

All  receivers  used  with 
noise,  15-40  dB  SNR, 

Monte  Carlo  simulation 
(600  iterations) 

trajallnmc.m 

SNR  dependent 

See  Figure  4.15 

SNR  dependent 
See  Figure  4.16 

Available  receivers  used 
with  no  noise 

tray group. m 

0.054 

2.200 

Available  receivers  used 
with  noise,  35  dB  SNR 

trajgroupn.m 

0.169 

7.681 

Available  receivers  used 
with  noise,  15-40  dB  SNR, 
Monte  Carlo  simulation 
(600  iterations) 

trajgroupnmc.  m 

SNR  dependent 

See  Figure  4.23 

SNR  dependent 
See  Figure  4.24 

Table  4.2.  Summary  of  TM  Tracker  Simulations 


53 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


54 


V.  CONCLUSION 


A.  SUMMARY 

The  objective  of  this  thesis  was  to  determine  system  requirements  and  limitations 
of  the  TM  Tracker  to  obtain  TSPI  with  1-meter  accuracy.  By  simulating  the  system 
under  ideal  conditions  with  no  noise  present,  a  theoretical  lower  limit  of  TSPI  accuracy 
was  established  using  all  remote  receivers  on  the  missile  range.  Noise  was  then  added  to 
Doppler  frequency  measurements  to  simulate  real-world  conditions  and  determine  effects 
on  TSPI  accuracy.  Based  on  Monte  Carlo  simulations,  if  all  nine  receivers  are  used 
throughout  the  entire  missile  trajectory  with  SNRs  above  approximately  39  dB,  the  TM 
Tracker  is  capable  of  producing  TSPI  accuracy  within  1  meter. 

Simulations  were  then  conducted  with  the  implementation  of  the  receivers’  45- 
degree,  3-dB  beam  widths.  Only  those  receivers  that  could  obtain  the  telemetry  signal 
within  their  respective  beam  widths  were  utilized  at  any  given  time  during  the  trajectory. 
Because  the  receiver  antennas  are  stationary,  only  one  window  of  time  could  be  used  for 
signal  reception  within  each  receiver’s  beam  width.  Special  consideration  had  to  be 
given  to  the  computation  of  the  initial  range  differences  when  a  receiver  was  added  after 
missile  launch.  Because  its  initial  range  difference  is  based  on  an  estimated  position,  to 
which  all  subsequent  differential  range  differences  are  added,  position  errors  after  the 
addition  of  a  receiver  during  the  trajectory  became  progressively  worse. 

The  receiver  addition  was  accomplished  late  in  the  missile’s  trajectory  to  prevent 
the  rapid  accumulation  of  errors.  Monte  Carlo  simulations  revealed  that  receivers 
required  SNRs  in  the  45-dB  range  to  obtain  1-meter  TSPI  accuracy.  Other  methods  of 
receiver  utilization  that  result  in  lower  errors  for  the  same  values  of  SNR  are  quite 
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possible,  but  an  enormous  number  of  combinations  exists,  and  each  Monte  Carlo 
simulation  takes  about  22  hours  to  complete  on  a  Pentium-Ill  class  computer. 

B.  FUTURE  RESEARCH 

A  zero-crossing  counter  approach  was  used  in  the  field-tested  TM  Tracker  for 
frequency  measurement,  which  was  then  converted  to  TDOA  and  a  range  difference. 
Sample  measurements  were  made  every  0.01  sec  and  consisted  of  the  number  of  zero 
crossings  from  time  start  to  the  measurement  point.  Because  the  count  started  at  time 
zero,  each  receiver  had  to  have  the  missile  within  full  view  during  its  entire  trajectory. 
MATLAB  simulations  of  the  TM  Tracker  assumed  frequency  measurements  were  limited 
by  the  theoretical  lower  bound  due  to  noise  and  did  not  take  into  account  the  implications 
of  using  the  zero-crossing  counter  approach  in  determining  frequency,  which  quantizes 
range  difference  measurements  to  the  order  of  a  wavelength,  (0.130  m).  [Ref.  7]  Because 
of  the  inherent  requirements  and  limitations  of  the  zero-crossing  counter  approach  and  its 
effects  on  TSPI  accuracy,  a  phase  detector  system  was  utilized  in  all  simulations.  The 
implementation  of  a  phase  detector  system  for  improved  frequency  measurement  and  a 
continuous  wave  (CW)  telemetry  signal  for  band  limiting  the  noise  within  the  remote 
receivers  is  currently  being  investigated  for  use  in  future  testing  of  the  TM  Tracker 
system. 
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APPENDIX  A.  MATLAB  SOURCE  CODES 


The  MATLAB  files  shown  in  Table  A.l  were  used  for  simulations  in  support  of 
this  thesis.  Their  respective  source  codes  are  listed  within  this  appendix. 


MATLAB  File 

Purpose 

fdoafang4nmc.  m 

Monte  Carlo  simulation  of  Fang’s  algorithm  for  determining  position 
error  vs.  SNR  for  simple  trajectory 

tdoafang4nmc.  m 

Monte  Carlo  simulation  of  Fang’s  algorithm  for  determining  X,  Y,  Z, 
and  total  position  error  vs.  RMS  measurement  error 

smith_abel.m 

Smith  -  Abel  algorithm  for  estimating  source  location 

trajall.m 

TM  Tracker  simulation  using  all  receivers  with  no  noise 

trajalln.m 

TM  Tracker  simulation  using  all  receivers  with  noise 

trajallnmc.m 

Monte  Carlo  simulation  of  TM  Tracker  using  all  receivers  with  noise 

angles.m 

Azimuth  and  elevation  angle  calculation  for  all  receivers 

traj group. m 

TM  Tracker  simulation  using  available  receivers  with  no  noise 

trajgroupn.m 

TM  Tracker  simulation  using  available  receivers  with  noise 

trajgroupnmc.  m 

Monte  Carlo  simulation  of  TM  Tracker  using  available  with  noise 

Table  A.l.  MATLAB  Files 
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fdoafang4nmc .  m 


%Position  from  FDOA  using  Fang  method 

%Four  receivers  at  ( 0 , 0 , 0 ) ,  (X2 , 0 , 0 ) ,  (X3 , Y3 , 0 ) ,  (X4 , Y4 , 0 ) 
%Inputs  source  position  and  calculates 
%Position  error  with  Monte  Carlo  simulation 
% f doaf ang4nmc . m 
clear,*clc;clf ; 

%Receiver  locations 


X2=-10 ;  %km 
X3=“10 ;  %km 
X4=-5 ;  %km 
Y3=5 ;  %km 
Y4=-5 ;  %km 


%Input  source  location 
xp=input('x  position^'); 
yp=input ( 'y  position= ' ) ; 
zp=input ( ' z  position= ' ) ; 

%Telemetry  and  missile  parameters 


f o=3e9 ;  %Hz 

v=341 . 6 ;  %m/s 

wl=3e8/fo;  %m 

f dm=v/wl ;  %Hz 

T=le-3 ;  %s 


%Generate  SNR  loop  (15*40  dB) 
for  snrdb=15 : 40 ; 

snr^lO^ (snrdb/10) ; 

%Theoretical  minimum  frequency  RMS  measurement  error 
dfrms=sqrt (3/snr) / (pi*T) ; 

%Generate  Monte  Carlo  loop  (10  iterations) 
for  mc=l : 10 ; 

%Compute  initial  ranges  from  receivers  to  source 
Rl=sqrt  (xp/s2+ypA2+zp/N2)  ; 

R2  =  sqrt  {  (xp-X2 )  A2+yp^2+zpA2  )  ; 

R3  =  sqrt  (  (xp-X3)  A2+ (yp-Y3 )  A2  +  zpA2 )  ; 

R4=sqrt  (  (xp-X4)  A2+ (yp*Y4)  A2+zpA2  )  ; 

%Compute  initial  RDs 

R2 1 i =R2 -R1 ; R2 1 =R2 1 i ; 

R3 li=R3 -R1 ; R3 1=R3 1 i ; 

R41i=R4-Rl ; R41=R41i ; 

%Define  initial  positions 
x=xp;xx=x; 

y=yp;yy=y; 
z=zp; zz=z ; 


%Define  time  base  for  trajectory 
for  t=0 : T : 10 ; 
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if  t==0;tx=0; 
else  tx=[tx  t];end; 

%Find  present  missile  position 
xpl=xp+ (v*t/1000) ; 
ypl=yp; 
zpl=zp; 

if  t==0;xpx=xpl;ypx=ypl;zpx=zpl; 

else  xpx=[xpx  xpl] ;ypx= [ypx  ypl ] ; zpx= [ zpx  zpl] ;end 

%Find  next  missile  position 
xm=xpl+v*T/1000 ; 
ym=ypl ; 
zm=zpl; 

%Find  Doppler  shifts  at  receivers 

%Present  missile  position 

phi0=atan2  (zpl,  sqrt  (xplA2+yplA2)  )  ; 
beta=atan2 (ypl , xpl ) ; 

phil=atan2 (zpl, sqrt ( (xpl-X2) A2+yplA2) ) ; 

alpha=atan2 (ypl, (xpl-X2) ) ; 

phi2=atan2 ( (ypl-Y3) , (xpl-X3) ) ; 

gamma=atan2 ( (ypl-Y4) , (xpl-X4) ) ; 

phi4=atan2 (zpl, sqrt ( (xpl-X3) A2+ (ypl-Y3 ) A2) ) ; 

phi5=atan2 (zpl, sqrt ( (xpl-X4) A2+ (ypl-Y4) A2) ) ; 

%Find  Doppler  frequencies  at  present  time 

fdl=fdm*cos (phiO) *cos (beta) +dfrms*randn; 
fd2=fdm*cos (phil) *cos (alpha) +dfrms*randn; 
fd3=fdm*cos (phi2 ) *cos (phi4) +dfrms*randn; 
fd4=fdm*cos (gamma) *cos (phi5) +dfrms*randn; 

%Find  differential  RDOAs  at  present  time 
drdoa3a=-le-4*T* (fdl-fd3)  ; 
drdoa2a=-le-4*T* (fdl-fd2) ; 
drdoa4a=-le-4*T* (fdl-fd4) ; 

%Find  Doppler  shifts  at  receivers 

%Next  missile  position 

phi0=atan2 (zm, sqrt (xmA2+ymA2) ) ; 
beta=atan2 (ym,xm) ; 

phil=atan2 (zm, sqrt( (xm-X2) A2+ymA2) ) ; 
alpha=atan2 (ym, (xm-X2) ) ; 
phi2=atan2 ( (ym-Y3) , (xm-X3) ) ; 
gamma=atan2 ( (ym-Y4) , (xm-X4) ) ; 
phi4=atan2 (zm, sqrt ( (xm-X3) A2+ (ym-Y3 ) A2 ) ) ; 
phi5=atan2 (zm, sqrt ( (xm-X4) A2+ (ym-Y4) A2 ) ) ; 

%Find  Doppler  frequencies  at  next  time 

fdl=fdm*cos (phiO) *cos (beta) +dfrms*randn; 
fd2=fdm*cos (phil) *cos (alpha) +dfrms*randn; 
fd3=fdm*cos (phi2 ) *cos (phi4) +dfrms*randn; 
fd4=fdm*cos (gamma) *cos (phi5) +dfrms*randn; 

%Find  differential  RDOAs  at  next  time 
drdoa3b=-le-4*T* (fdl-fd3) ; 
drdoa2b=-le-4*T* (fdl-fd2) ; 
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drdoa4b=-le-4*T* (fdl-fd4)  ; 


%Find  average  differential  RDOAs 

drdoa2= (drdoa2a+drdoa2b) /2; 
drdoa3= (drdoa3a+drdoa3b) /2; 
drdoa4= (drdoa4a+drdoa4b) /2 ; 

%Find  RDOAs 

R21a=R21+drdoa2 ; 

R31a=R31+drdoa3 ; 

R41a=R41+drdoa4 ; 

%Find  coefficients 

gl= (R31a*X2/R21a-X3) /Y3 ; 
g2= (R41a*X2/R21a-X4) /Y4; 
g=gl-g2; 

hl=  (X3A2+Y3 A2-R31aA2+R31a*R21a*  (1-  (X2/R21a)  A2) )  /  (2*Y3) 
h2= (X4A2+Y4A2-R41aA2+R41a*R21a* (1- (X2/R21a) A2) ) / (2*Y4) 
h=hl-h2 ; 

dl=- ( (1-(X2/R21a) A2)+gl*gl) ; 
d2=- ( (1-(X2/R21a) A2)+g2*g2)  ; 
d=dl-d2 ; 

el=X2*( (l-(X2/R21a)A2) )-2*gl*hl; 
e2=X2*( ( 1- (X2/R21a) A2 ) )-2*g2*h2; 
e=el-e2 ; 

f 1= ( ( R21aA2 ) /4 ) * ( 1- (X2 /R21a ) A2 ) A2-hl*hl ; 
f2= ( (R21aA2 ) /4) * ( 1- (X2/R21a) A2 ) A2-h2*h2 ; 
f =f 1-f 2 ; 

%Reset  range  differences 
R21=R21a; 

R31=R31a; 

R41=R41a; 

%Find  new  positions 

x= (-e+sqrt (e*e-4*d*f ) ) / (2*d) ; 
y=gl*x+hl ; 

za=sqrt(dl*x*x+el*x+fl)  ; 
z=abs (za) ; 

%Find  position  error  in  meters 

err=1000*sqrt ( (x-xm) A2+ (y-ym) A2+ (z-zm) A2) ; 

%Fill  error  matrix  for  missile  trajectory 
if  t==0;xe=err; 
else  xe=[xe  err]; end 

% Re turn  to  trajectory  loop 
end; 

%Calculate  mean  error  for  missile  trajectory 
merr=mean (xe) ; 

%Fill  mean  error  matrix  for  Monte  Carlo  simulation 
if  mc==l ;mcxe=merr ; 
else  mcxe=[mcxe  merr] ;end 
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%Return  to  Monte  Carlo  loop 
end; 

%Calculate  mean  error  for  Monte  Carlo  iterations 
mcmerr==mean  (mcxe)  ; 

%Fill  mean  error  matrix  for  each  value  of  SNR 
if  snrdb==15 ; snmcmerr=mcmerr ; snrdbm=snrdb; 
else  snmcmerr= [snmcmerr  mcmerr] ; snrdbm= [snrdbm  snrdb];end 

%Output  results 

fprintf ( '  SNR  =  %g\n'  ,  snrdb) ; 
fprintf ( '  x  =  %g  \n'  ,  x)  ; 

fprintf ( '  y  =  %g  \n'  ,  y)  ; 

fprintf ( '  z  =  %g  \n'  ,  z) ; 

fprintf ( '  Mean  error  for  10  iterations  =  %g  \n'  ,  mcmerr) ; 

%Return  to  SNR  loop 
end; 

%Plot  the  graph 
figured); 

plot (snrdbm, snmcmerr) ;grid; 
title (' Position  Error  vs.  SNR'); 

xlabel ( ' SNR  -  dB ') ;ylabel (' Position  Error  -  meters'); 


t  doafang4nmc .  an 

%Position  from  TDOA  using  Fang  method 

%Four  receivers  at  (0, 0, 0) , (X2, 0, 0) , (X3, Y3, 0) , (X4, Y4, 0) 

% Inputs  source  position,  adds  TDOA  measurement  errors,  and  calculates 
%Position  error  with  Monte  Carlo  simulation 
%tdoafang4nmc .m 
clear;clc;clf ; 

%Receiver  locations 


X2=-10;  %km 
X3=-10 ;  %km 
X4=-5 ;  %km 
Y3=5 ;  %km 
Y4=-5 ; 


% Input  source  location 
xp=input('x  position='); 
yp=input('y  position='); 
zp=input('z  position='); 

%Compute  initial  ranges  from  receivers  to  source 
Rl=sqrt  (xp/^2+yp/s2+zp^2 )  ; 

R2=sqrt  (  (xp-X2)  /N2+yp^2  +  zp^2 )  ; 

R3=sqrt ( (xp-X3) ^2+ (yp-Y3) A2+zp^2) ; 

R4=sqrt ( (xp-X4 ) A2+ (yp-Y4) A2+zpA2) ; 

%Generate  TDOA  error  loop 
for  j=. 000001: .0000005: .00001; 
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%Generate  Monte  Carlo  loop  (10,000  iterations) 
for  i=l: 10000; 

%Compute  TDOA 

R21= (R2-R1 ) * (1+ j  *randn) ; 

R31= (R3-R1 ) * ( 1+ j  *randn) ; 

R41= (R4-R1 ) * (1+ j  *randn) ; 

%Find  coefficients 

gl=(R31*X2/R21-X3) /Y3; 
g2= (R41*X2/R21-X4) /Y4; 
g=gl-g2; 

hl=  (X3 A2+Y3A2-R31A2+R31*R21*  (1-  (X2/R21)  A2)  )  /  (2*Y3)  ; 
h2= (X4 A2+Y4 A2-R41 A2+R41*R21* (1-(X2/R21) A2) ) / (2*Y4)  ; 
h=hl-h2 ; 

dl=- ( (1- (X2/R21) A2) +gl*gl) ; 
d2=- ( (1- (X2/R21) "2) +g2*g2) ; 
d=dl-d2 ; 

el=X2* ( (1- (X2/R21) A2) ) -2*gl*hl; 
e2=X2* ( (1- (X2/R21) A2) ) -2*g2*h2; 
e=el-e2 ; 

fl= ( (R21A2) /4) * (1- (X2/R21) A2) A2-hl*hl; 
f 2= ( (R21A2) /4) * (1-(X2/R21) A2) A2-h2*h2; 
f =f 1-f 2 ; 

%Find  positions 

x(i) = (-e+sqrt (e*e-4*d*f ) ) / (2*d) ; 
y (i) =gl*x(i) +hl; 

z (i) =sqrt (dl*x(i) *x(i) +el*x(i) +fl)  ; 

%Find  position  errors  in  meters 

xerr(i)  =  1000*abs (xp-x(i) ) ; 
yerr(i)  =  1000*abs (yp-y (i) ) ; 
zerr(i)  =  1000*abs (zp-z (i) ) ; 

err(i) =1000*sqrt( (x(i) -xp) A2+(y(i) -yp) A2+(z (i) -zp) A 

%Retum  to  Monte  Carlo  loop 
end; 

%Find  measurement  errors  in  inches 
R21p= (R2-R1) * j *1000*39. 37; 

R31p= (R3-Rl)*j *1000*39. 37; 

R41p= (R4-R1) * j*1000*39 .37; 

%Find  and  calculate  smallest  RMS  measurement  error 
if  R21p<R31p  &  R21p<R41p; 

RMS=sqrt ( (R21p) A2) ; 
elseif  R31p<R21p  &  R31p<R41p; 

RMS=sqrt ( (R31p) A2) ; 
else  RMS=sqrt ( (R41p) A2 )  ; 
end; 


%0utput  results 

fprintf (' Introduced  TDOA  Error  (km) 
fprintf('RMS  TDOA  Measurement  Error 
fprintf ( 'Mean  x-error  -  meters  =  %g 
fprintf ( 'Mean  y-error  -  meters  =  %g 


=  %g\n'  ,  j) ; 

-  inches  =  %g  \n'  , 
\n'  ,  mean(xerr)); 
\n'  ,  mean(yerr)); 


RMS) 
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fprintf ( 'Mean  z-error  -  meters  =  %g  \n'  ,  mean(zerr) ) ; 
fprintf ( 'Mean  total  error  -  meters  =  %g  \n'  ,  mean(err)); 
fprintf ( '  \n' )  ; 

%Fill  matrices  for  plotting  values 
if  j  == . 000001 ; 

nn=RMS;xx=mean(xerr) ;yy=mean (yerr) ; zz=mean ( zerr) ;ee=mean(err) ; 
else  nn= [ nn  RMS];xx=[xx  mean (xerr) ] ;yy= [yy  mean(yerr)]; 

zz=[zz  mean ( zerr) ]; ee= [ee  mean (err) ]; 
end; 

%Return  to  TDOA  error  loop 
end; 

%Plot  the  graphs 
figured)  ; 
plot(nn,xx) ;grid; 

title ( 'X-Position  Error  vs.  RMS  Measurement  Error'); 
xlabel ( 'RMS  Measurement  Error  -  inches'); 
ylabel ( 'X-Position  Error  -  meters'); 

figure (2 ) ; 
plot (nn,yy) ;grid; 

title ( 'Y-Position  Error  vs.  RMS  Measurement  Error'); 
xlabel ('RMS  Measurement  Error  -  inches'); 
ylabel ( 'Y-Position  Error  -  meters'); 

figured)  ; 

plot (nn, zz) ; grid; 

title ( 'Z-Position  Error  vs.  RMS  Measurement  Error'); 
xlabel ( ' RMS  Measurement  Error  -  inches ' ) ; 
ylabel ( 'Z-Position  Error  -  meters'); 

figure (4) ; 

plot (nn, ee) ;grid; 

title( 'Total  Position  Error  vs.  RMS  Measurement  Error'); 
xlabel ( ' RMS  Measurement  Error  -  inches ' ) ; 
ylabel ( 'Total  Position  Error  -  meters'); 


smith_abel  .m 

%Smith-Abel  algorithm  used  as  a  function 

%Computes  source  location  estimate 

%Inputs  are  receiver  locations,  range  differences, 

%Weights,  and  receiver  ranges  to  reference  receiver 
%Algorithm  is  tailored  for  nine  receivers 
%smith_abel .m 

function  [source_estimate] =smith_abel (RL, d, w,  R) 

%Source  locations  of  receivers  other  than  reference  receiver 
S=RL (2:9, : ) ; 

%Translate  coordinate  system  so  reference  receiver  is  at  origin 
S ( : , 1) =S ( : , 1 ) -RL (1,1) ;  %x-dimension  translated 
S ( : , 2 ) =S ( : , 2 ) -RL (1,2);  %y-dimension  translated 
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S(:,3)=S(:,3)-RL(1,3)  ; 


%z-dimension  translated 


%Define  delta  matrix 
del ta=R . A2 -d . A2 ; 

%Define  8x8  diagonal  matrix 
W=diag (w) ; 

%Utilize  spherical  interpolation 
%Define  Ps  and  Pso  matrices 
Ps  =  S*inv(S'*W*S) *S'*W; 

Pso=eye(8)  -  Ps; 

%Compute  minimizing  value  for  Rs 

Rshat= (d'*Pso*W*Pso*delta) / {2*d' *Pso*W*Pso*d) ; 

%Compute  source  location  estimate 

estimate= . 5*inv ( S ' *W *S)  *S ' *W* (delta-2 *Rshat*d) ; 

%Translate  back  to  orignial  coordinate  system 
source__estimate=estimate '  +RL  ( 1 ,  :  )  ; 


trajall  .jn 

%Position  from  FDOA  using  Smith-Abel  algorithm 
%Inputs  missile  trajectory  and  receiver  SNR 

%Calculates  RDOA  from  FDOA  and  position  with  SNR  measurement  error 

%A11  receivers  utilized  throughout  entire  trajectory 

%trajall.m 

%Load  White  Sands  missile  trajectory 
load  traj ; 

t=traj  (1,0  ;xm=traj  (2,  : )  ;ym=traj  (3,  : )  ;  zm=tra  j  (4,  : )  ; 


% Smooth  missile  trajectory 
xms=polyval  (polyfit  (t  ,xm,  14)  ,  t)  ; 
yms=polyval  (polyfit  (t  ,ym,  14)  ,  t)  ; 
zms=polyval  (polyf  it  ( t ,  zm,  14)  ,  t )  ; 

%Obtain  polynomial  coefficients  for  trajectory 
xpc=polyf it ( t/Xm, 14)  ; 
ypc=polyfit ( t/ym, 14)  ; 
zpc=polyf it { t ,  zm,  14)  ; 

%Def ine  coefficients  for  velocity  (first  derivative  of  polynomials) 
i=l : 14 ; 

xvc= ( 15-i) . *xpc (i) ; 
yvc= (15-i) . *ypc (i) ; 
zvc= (15-i) . *zpc (i) ; 

% Evaluate  velocity  polynomials  to  obtain  velocity  component  vectors 
vxs=polyval (xvc, t) ; 
vys=polyval (yvc7 1) ; 
vzs=polyval  (zvc,  t)  ; 
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%White  Sands  remote  receiver  positions  (feet) 

RLi=  [ -  9243 .4  -3085.7  -127;  -9193.6  -3072.3  -125;  -9278.7  -3201.8 

-92.1;  -8992.9  -2879.8  -109.3;  1757  430.4  40.4;  -6458.5  -2112.3 
-39.3;  -3704.2  -1214.7  -33.6;  607  -3517.8  208.3;  -5305.4  2653 
-153.6] ; 


%Initialize  range  difference  vectors 

d=zeros (8,1) ; 

da=zeros (8,1)  ; 

di=zeros (8,1) ; 

db=zeros (8,1) ; 

% Range  from  non-reference  receiver  to  reference  receiver 
R=zeros (8,1) ; 


R(l) 

= 

norm(RL(l, : 

)  -  RL (2  , 

:)) 

R(2) 

= 

norm(RL(l, : 

)  -  RL (3  , 

:)  ) 

R(3) 

= 

norm(RL(l, : 

)  -  RL  (4 , 

:)  ) 

R(4) 

norm(RL(l, : 

)  -  RL ( 5 / 

:)  ) 

R{5) 

= 

norm(RL(l, : 

)  -  RL  ( 6 , 

:)) 

R<6) 

= 

norm(RL (1 , : 

)  -  RL  (7 , 

:)) 

R(7) 

= 

norm  (RL ( 1 ,  : 

)  -  RL  ( 8 , 

:)  ) 

R  (8) 

= 

norm(RL(l, ; 

)  -  RL (9, 

:)) 

%Set  all  weights  equal  to  one  (using  all  receivers) 
w=ones (8,1) ; 

%Telemetry  and  missile  data 


f o=2 . 3e9 ;  %  Hz 

wl= (3e8/fo) /.3048;  %  feet 

T=le-2 ;  %  sec 

Tmax=6.41;  %  sec 


% Initial  range  differences,  with  Receiver  1  as  the  reference 
sourcei= [xms (1)  yms(l)  zms(l)]; 


did) 
di  (2) 
di  (3 ) 
di  (4 ) 
di  (5) 
di  ( 6 ) 
di  (7 ) 
di  ( 8 ) 


norm(sourcei 
norm(sourcei 
norm(sourcei 
norm(sourcei 
norm (sourcei 
norm(sourcei 
norm  (sourcei 
norm (sourcei 


RL  (2 , 

:))  - 

RL  (3 , 

:))  - 

RL  (4, 

:))  - 

RL(5, 

:))  - 

RL  ( 6 , 

:))  - 

RL  (7 , 

:  )  )  - 

RL  ( 8 , 

:))  - 

RL  ( 9 , 

:))  - 

norm (sourcei 
norm (sourcei 
norm  (sourcei 
norm  (sourcei 
norm (sourcei 
norm (sourcei 
norm (sourcei 
norm (sourcei 


-  RL ( 1 ,  :  )  ) 

-  RL ( 1 ,  :)) 

-  RL ( 1 ,  :  )  ) 

-  RL ( 1 ,  :)  ) 

-  RL ( 1 , : ) ) 

-  RL ( 1 , : ) ) 

-  RL ( 1 , : ) ) 

-  RL ( 1 ,  :  )  ) 


d=di ; 


% Input  receiver  SNR 
snrdb= input ( ' snr (dB)  = ' )  ; 
snr=l(d  (snrdb/10) ; 

%Theoretical  minimum  frequency  RMS  measurement  error 
dfrms=sqrt (3/snr) / (pi*T) ; 


%Generate  loop  for  sampling  trajectory  points 

%Read  positions  and  velocities  at  each  sample  point  (every  .01  sec) 
for  ts=0 : T : Tmax; 
sample=ts/T+l; 
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sample=round (sample) ;  %Sample  must  be  an  integer  value 

x=xms (1, sample) ;y=yms (1, sample) ; z=zms  (1, sample) ; 

xl=xms ( 1 , sample+1 ) ;yl=yms ( 1 , sample+1) ; zl=zms ( 1 , sarrtple+l ) ; 

v= [vxs (1, sample)  vys (1, sample)  vzs (1, sample) ] ; 

vl= [vxs ( 1 , sample+1 )  vys ( 1 , sample+1 )  vzs ( 1 , sample+1 ) ] ; 

%Find  source  position  at  present  time 
source= [x  y  z]  ; 

%Find  source  position  at  next  time 
sourcel=[xl  yl  zl] ; 

%Find  source-to-receiver  vectors  at  present  time 
Bl=source-RL ( 1 , : ) ; 

B2=source-RL (2 , : ) ; 

B3=source-RL (3 , : ) ; 

B4=source-RL(4, : ) ; 

B5=source-RL (5 , : ) ; 

B6=source-RL ( 6 , : ) ; 

B7=source-RL (7 , : ) ; 

B8=source-RL ( 8 , : ) ; 

B9=source-RL ( 9 , : ) ; 

%Find  Doppler  frequency  vector  at  present  time 
F=v/wl ; 

%Find  Doppler  frequencies  at  receivers  at  present  time 
fdl=dot (Blf F) / norm ( B1 ) +df rms ; 
fd2=dot (B2 , F) /norm(B2 ) +dfrms ; 
f d3=dot (B3 , F) /norm(B3)+dfrms; 
fd4=dot (B4 , F) /norm(B4) +dfrms; 
fd5=dot (B5 , F) /norm(B5 ) +df rms ; 
fd6=dot (B6 , F) /norm(B6)+dfrms; 
fd7=dot  (B7,F)  /no2rm(B7  )  +dfrms; 
f d8=dot ( B8 , F) /norm (B8 ) +df rms; 
fd9=dot (B9 , F) /norm(B9) +dfrms; 

%Find  differential  RDOAs  at  present  time 
drdoa2a=-wl*T* (fdl-fd2) ; 
drdoa3a=-wl*T* (fdl-fd3 ) ; 
drdoa4a=-wl*T* (fdl-fd4) ; 
drdoa5a=-wl*T* (fdl-fd5) ; 
drdoa6a=-wl*T* (fdl-fd6) ; 
drdoa7a=-wl*T* (fdl-fd7 ) ; 
drdoa8a=-wl,lrT*  (fdl-fd8)  ; 
drdoa9a=-wl*T* (fdl-fd9) ; 

%Find  source-to-receiver  vectors  at  next  time 
Bl=sourcel-RL(l,  : )  ; 

B2=sourcel-RL (2 ,:); 

B3=sourcel-RL (3 , : ) ; 

B4=sourcel-RL (4 ,  : )  ; 

B5=sourcel-RL(5,  : )  ; 

B6  =  sourcel-RL ( 6 ,  : )  ; 

B7=sourcel-RL (7,  : )  ; 

B8  =  sourcel-RL(8/  ; )  ; 

B9=sourcel-RL (9, : ) ; 
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%Find  Doppler  frequency  vector  at  next  time 
Fl=vl/wl? 

%Find  Doppler  frequencies  at  receivers  at  next  time 
fdl=dot (Bl, FI) /norm (B1 ) +df rms ; 
fd2=dot (B2,Fl) /norm (B2 ) +df  rms ; 
fd3=dot (B3 , FI) /norm ( B3 ) +df rms; 
fd4=dot  (B4, FI)  /norm (B4 )  +df rms; 
fd5=dot (B5, FI) /norm (B5 ) +df  rms; 
fd6=dot (B6 , FI) /norm{B6) +dfrms; 
fd7=dot (B7 , FI ) /norm (B7 ) +df rms ; 
f d8=dot ( B8 , FI ) /norm ( B8 ) +df  rms ; 
fd9=dot ( B9 , FI ) /norm ( B9 ) +df  rms ; 

%Find  differential  RDOAs  at  next  time 
drdoa2b=-wl*T* (fdl-fd2) ; 
drdoa3b=-wl*T* (fdl-fd3) ; 
drdoa4b=-wl*T* (fdl-fd4) ; 
drdoa5b=-wl*T* (fdl-fd5) ; 
drdoa6b=-wl*T* (fdl-fd6)  ; 
drdoa7b=-wl *T* (fdl-fd7) ; 
drdoa8b=-wl*T* (fdl-fd8) ; 
drdoa9b=-wl*T* (fdl-fd9) ; 


%Find  average  differential  RDOAs 
drdoa2= (drdoa2a+drdoa2b) /2; 
drdoa3= (drdoa3a+drdoa3b) /2 ; 
drdoa4= (drdoa4a+drdoa4b) 12; 
drdoa5= (drdoa5a+drdoa5b) /2 ; 
drdoa6= (drdoa6a+drdoa6b) 12; 
drdoa7= (drdoa7a+drdoa7b) /2; 
drdoa8= (drdoa8a+drdoa8b) 12; 
drdoa9= (drdoa9a+drdoa9b) / 2; 


%Find  RDOAs 

da ( 1 ) =d ( 1 ) +drdoa2 ; 
da ( 2 ) =d { 2 ) +drdoa3 ; 
da ( 3 ) =d ( 3 ) +drdoa4 ; 
da ( 4 ) =d { 4 ) +drdoa5 ; 
da ( 5 ) =d { 5 ) +drdoa6 ; 
da ( 6 ) =d ( 6 ) +drdoa7 ; 
da ( 7 ) =d ( 7 ) +drdoa8 ; 
da ( 8 ) =d ( 8 ) +drdoa9 ; 

%Reset  range  differences 
d=da; 


%Find  RDOAs  from  geometry 


db(l) 

=  normfsourcel 

-  RL  (2 , 

:)  ) 

db(2) 

=  normfsourcel 

-  RL ( 3 , 

:)  ) 

db  (3 ) 

=  normfsourcel 

-  RL (4 , 

:)  ) 

db(4) 

=  normfsourcel 

-  RL ( 5 , 

:)  ) 

db(5) 

=  normfsourcel 

-  RL ( 6 , 

:)  ) 

db(6) 

=  normfsourcel 

-  RL ( 7 , 

:  )  ) 

db(7) 

=  normfsourcel 

-  RL(8, 

:  )  ) 

db(8) 

=  normfsourcel 

-  RL ( 9 , 

:)  ) 

norm(sourcel  - 
norm(sourcel  - 
norm(sourcel  - 
normfsourcel  - 
norm(sourcel  - 
normfsourcel  - 
normfsourcel  - 
normfsourcel  - 


RL(1, :) ) 
RL ( 1 , : ) ) 
RL(1, : ) ) 
RLfl, :) ) 
RLfl, :) ) 
RLfl, :) ) 
RLfl, :) ) 
RLfl, : ) ) 
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%Fill  matrices  for  time,  position,  and  velocity  components 
if  ts==0;  t=ts+T;xx=x;yy=y;  zz=z;vx=v(l)  ;vy=v(2)  ;vz=v(3)  ; 
else  xx=  [xx  x]  ; yy=  [yy  y] .;  zz  =  [zz  z]  ;  t  =  [t  ts+T] ; . . . 
vx=[vx  v ( 1 ) ] ; vy= [vy  v(2) ] ;vz=[vz  v(3) ] ; 

end; 

%Estimated  RDOA  error 
errt=norm(d  -  db)  ; 

%Fill  matrix  for  estimated  RDOA  error 
if  ts==0 ;xet=errt ; 
else  xet=[xet  errt] ; 
end; 

%Apply  Smith-Abel  algorithm  to  find  estimated  position  from  RDOAs 
source_est=smith__abel  (RL,  d,  w,  R)  ; 

%Estimated  position  error 

errp=norm(sourcel  -  source_est) ; 

%Fill  matrix  for  estimated  position  error 
if  ts==0 ; xep=errp ; 
else  xep=[xep  errp] ; 
end 

%Return  to  trajectory  loop 
end ; 

%Obtain  position  and  velocity  components  vectors 
pos=[t'  xx'  yy'  zz']; 
pos=pos ' ; 

v= [ t '  vx '  vy '  vz ' ] ; 
v=v '  ; 

%Output  results  to  screen 
merrt=mean (xet ) ; 

fprintf  ( 'Mean  RDOA  Error  (ft)  =  %g  \n'  ,merrt); 
merrp=mean(xep) ; 

fprintf ( 'Mean  Position  Error  (ft)  =  %g  \n'  ,merrp) ; 

%Plot  graphs 

figure (1) 

plot (t, xet) ; grid; 

title ('RDOA  Error  vs  Time'); 

xlabel('Time  -  sec ') ;y label ( 'RDOA  error  -  ft'); 

figure (2) 

plot(t,xep) ; grid; 

title (' Position  Error  vs  Time'); 

xlabel('Time  -  sec ')  ,-ylabel  ('  Position  Error  -  ft'); 
figure (3) 

plot 3 (xx,yy, zz) ; grid; 
hold  on 

plot3  (RL(  :  ,1)  ,RL(  :  ,2)  ,  RL  (  :  ,3)  ,  '*')  ; 
hold  off 
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title ( 'Missile  Trajectory  and  Remote  Receiver  Locations'); 
xlabel ( ' f t ' ) ;y label ( ' f t ' ) ; zlabel ( ' f t ' ) ; 


trajalln.m 

%Position  from  FDOA  using  Smith-Abel  algorithm 

%Inputs  missile  trajectory  and  receiver  SNR  and  calculates 

%RDOA  from  FDOA  and  position  with  SNR  measurement  error 

%Random  white  Gaussian  noise  added  to  signal 

%A11  receivers  utilized  throughout  entire  trajectory 

%trajalln.m 

%Load  White  Sands  missile  trajectory 
load  tra j ; 

t=traj (1, : ) ;xm=traj (2, : ) ;ym=traj (3, : ) ;zm=traj (4, : ) ; 

%Smooth  missile  trajectory 
xms=polyval (polyf it (t ,  xm, 14)  ,  t) ; 
yms=polyval (polyf it (t , ym, 14) , t ) ; 
zms=polyval (polyf it (t, zm, 14) , t ) ; 

%Obtain  polynomial  coefficients  for  trajectory 
xpc=polyfit (t ,xm, 14) ; 
ypc=polyf it (t,ym, 14) ; 
zpc=polyf it (t, zm, 14) ; 

%Def ine  coefficients  for  velocity  (first  derivative  of  polynomials) 
i=l : 14 ; 

xvc= (15-i) . *xpc (i)  ; 
yvc=  (15-i)  .  *ypc  (i)  ; 
zvc= (15-i) . *zpc (i) ; 

% Evaluate  velocity  polynomials  to  obtain  velocity  component  vectors 
vxs=polyval (xvc , t ) ; 
vys=polyval (yvc, t) ; 
vzs=polyval ( zvc  #  t ) ; 

% White  Sands  remote  receiver  positions  (feet) 

RL= [ -9243 . 4  -3085.7  -127;  -9193.6  -3072.3  -125;  -9278.7  -3201.8 

-92.1;  -8992.9  -2879.8  -109.3;  1757  430.4  40.4;  -6458.5  -2112.3 
-39.3;  -3704.2  -1214.7  -33.6;  607  -3517.8  208.3;  -5305.4  2653 
-153.6] ; 

%Initialize  range  difference  vectors 

d=zeros (8,1) ; 

da=zeros (8,1) ; 

di=zeros (8,1) ; 

db=zeros (8,1) ; 

%Range  from  non-reference  receiver  to  reference  receiver 
R=zeros (8,1) ; 

R ( 1 )  =  norm (RL ( 1 , : )  -  RL(2,:)); 

R ( 2 )  =  norm(RL(l,:)  -  RL{3,:)); 

R (3 )  =  norm(RL(l, : )  -RL(4,:)); 

R (4 )  =  norm (RL ( 1 , : )  -  RL(5,:)); 


R ( 5 )  =  norm(RL (1, : )  -  RL(6,:)); 

R ( 6 )  =  norm(RL(l, :)  -RL(7,:)); 

R ( 7 )  =  norm(RL(l, :)  -RL(8,:)); 

R ( 8 )  =  norm(RL(l, : )  -  RL(9,:)); 

%Set  all  weights  equal  to  one  (using  all  receivers) 
w=ones (8,1); 

%Telemetry  and  sampling  data 


fo=2 . 3e9 ;  %  Hz 

wl=(3e8/fo)/.3048;  %  feet 

T=le-2;  %  sec 

Tmax=6.41;  %  sec 


%Initial  range  differences,  with  Receiver  1  as  the  reference 
sourcei= [xms (1)  yms(l)  zms(l)]; 

di(l)  =  norm(sourcei  -  RL(2,:))  -  norm(sourcei  -  RL(1,:)); 

di (2)  =  norm(sourcei  -  RL(3,:))  -  norm(sourcei  -  RL(1,:)); 

di (3 )  =  norm(sourcei  -  RL(4,:))  -  norm(sourcei  -  RL(1,:)); 

di (4)  =  norm(sourcei  -  RL(5,:))  -  normfsourcei  -  RL(1,:)); 

di (5)  =  norm(sourcei  -  RL(6,:))  -  norm(sourcei  -  RL(1,:)); 

di(6)  =  norm(sourcei  -  RL(7,:))  -  norm(sourcei  -  RL(1,:)); 

di (7)  =  norm( sourcei  -  RL(8,:))  -  norm(sourcei  -  RL(1,:)); 

di ( 8 )  =  norm (sourcei  -  RL(9,:))  -  norm (sourcei  -  RL(1,:)); 

d=di ; 

% Input  receiver  SNR 
snrdb= input ( ' snr (dB) = ' )  ; 
snr=10A (snrdb/10) ; 

%Theoretical  minimum  frequency  RMS  measurement  error 
dfrms=sqrt (3/snr) / (pi*T) ; 

%Generate  loop  for  sampling  trajectory  points 

%Read  positions  and  velocities  at  each  sample  point  (every  .01  sec) 
for  ts=0 : T : Tmax; 
sample=ts/T+l; 

sample=round ( sample) ;  %Sample  must  be  an  integer  value 

x=xms (1, sample) ;y=yms  (1, sample) ; z=zms (1, sample) ; 

xl=xms (1, sample+1) ;yl=yms ( 1 , sample+1) ; zl=zms (1, sample+1) ; 

v= [vxs (1, sample)  vys (1, sample)  vzs (1/ sample) ] ; 

vl= [vxs (1, sample+1)  vys (1, sample+1)  vzs (1, sample+1) ] ; 

%Find  source  position  at  present  time 
source=[x  y  z] ; 

%Find  source  position  at  next  time 
sourcel=[xl  yl  zl] ; 

%Find  source-to-receiver  vectors  at  present  time 
Bl=SOurce-RL ( 1 , : ) ; 

B2=source-RL (2 , : ) ; 

B3 = source -RL (3 , :) ; 

B4=source-RL(4, : ) ; 

B5=source-RL(5, :) ; 
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time 


B6=source-RL ( 6 ,  :  )  ; 

B7=source-RL (7,  :  )  ; 

B8=source-RL ( 8  ,  : )  ; 

B9=source-RL(9,  : )  ; 

%Find  Doppler  frequency  vector  at  present 
F=v/wl ; 

%Find  Doppler  frequencies  at  receivers  at  present  time  (with  noise) 
fdl^dot ( B1 , F) /norm (B1 ) +df rms*randn; 
fd2=dot (B2,  F) /norm(B2) +dfrms*randn; 
fd3=dot (B3 , F) /norm(B3) +dfrms*randn; 
fd4=dot (B4, F) /norm(B4) +dfrms*randn; 
fd5=dot (B5 , F) /norm (B5 ) +dfrms*randn; 
fd6=dot (B6 , F) /norm (B6 ) +dfrms*randn; 
fd7=dot (B7 , F) /norm (B7 ) +dfrms*randn; 
f d8=dot (B8 , F) /norm (B8 ) +df rms*randn; 
fd9=dot (B9 , F) /norm(B9) +dfrms*randn; 

%Find  differential  RDOAs  at  present  time 
drdoa2a=~wl*T* (fdl-fd2)  ; 
drdoa3a=-wl*T* (fdl-fd3)  ; 
drdoa4a=-wl*T* (fdl-fd4)  ; 
drdoa5a=-wl*T* (fdl-fd5)  ; 
drdoa6a=-wl*T* (fdl-fd6)  ; 
drdoa7a=-wl*T* (fdl-fd7)  ; 
drdoa8a=-wl*T* (fdl-fd8)  ; 
drdoa9a=-wl*T* (fdl-fd9)  ; 

%Find  source-to-receiver  vectors  at  next  time 
Bl=sourcel-RL(l, : ) ; 

B2=sourcel-RL (2 , : )  ; 

B3=sourcel-RL (3 , :  )  ; 

B4=sourcel-RL(4,  :  )  ; 

B5=sourcel-RL(5,  :  )  ; 

B6=sourcel-RL ( 6 ,  :  )  ; 

B7=sourcel-RL (7  ,  :  )  ? 

B8=sourcel-RL(8,  :  )  ; 

B9=sourcel-RL (9 ,  :  )  ? 

%Find  Doppler  frequency  vector  at  next  time 
Fl=vl/wl ; 

%Find  Doppler  frequencies  at  receivers  at  next  time  (with  noise) 
fdl=dot (Bl, FI) /norm(Bl) +dfrms*randn; 
f d2=dot (B2 , FI) /norm(B2 ) +df rms*randn; 
fd3=dot (B3 , FI ) /norm(B3 ) +df rms*randn; 
f d4=dot (B4, FI) /norm(B4) +dfrms*randn; 
f d5=dot (B5 , FI ) /norm(B5 ) +df rms*randn; 
f d6=dot (B6,F1) /norm(B6) +dfrms*randn; 
fd7=dot (B7 , FI) /norm(B7) +dfrms*randn; 
fd8=dot (B8,F1) /norm(B8)+dfrms*randn; 
fd9=dot (B9, FI) /norm(B9) +dfrms*randn; 

%Find  differential  RDOAs  at  next  time 
drdoa2b=-wl*T* ( fdl-f d2  )  ; 
drdoa3b=-wl*T* ( fdl-fd3 )  ; 
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drdoa4b=-wl*T* (fdl-fd4)  ; 
drdoa5b=-wl*T* (fdl-fd5)  ; 
drdoa6b=-wl*T* (fdl-fd6)  ; 
drdoa7b=-wl*T* (fdl-fd7)  ; 
drdoa8b=-wl*T* (fdl-fd8)  ; 
drdoa9b=-wl*T* (fdl-fd9)  ; 


%Find  average  differential  RDOAs 
drdoa2= (drdoa2a+drdoa2b) /2 ; 
drdoa3= (drdoa3a+drdoa3b) / 2 ; 
drdoa4= (drdoa4a+drdoa4b) /2 ; 
drdoa5= (drdoa5a+drdoa5b) / 2 ; 
drdoa6= (drdoa6a+drdoa6b) / 2 ; 
drdoa7= (drdoa7a+drdoa7b) /2 ; 
drdoa8= (drdoa8a+drdoa8b) / 2 ; 
drdoa9= (drdoa9a+drdoa9b) / 2 ; 


%Find  RDOAs 

da ( 1 ) =d ( 1 ) +drdoa2 ; 
da ( 2 ) =d ( 2 ) +drdoa3 ; 
da ( 3 ) =d ( 3 ) +drdoa4 ; 
da ( 4 ) =d { 4 ) +drdoa5 ; 
da ( 5 ) =d ( 5 ) +drdoa6 ; 
da ( 6 )  =d ( 6 ) +drdoa7 ; 
da ( 7 ) =d ( 7 ) +drdoa8 ; 
da ( 8 ) =d ( 8 ) +drdoa9 ; 

%Reset  range  differences 
d=da; 


%Find  RDOAs  from  geometry 
db(l)  =  norm(sourcel  - 
db(2)  =  norm(sourcel  - 

db(3)  =  norm(sourcel  - 
db(4)  =  norm(sourcel  - 
db(5)  =  norm(sourcel  - 

db ( 6 )  =  norm(sourcel  - 
db(7)  =  norm(sourcel  - 
db(8)  =  norm(sourcel  - 
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:)); 

RL  ( 1 , 

:)); 

RL  ( 1 , 

:)); 

%Fill  matrices  for  time,  position,  and  velocity  components 
if  ts==0;t=tS+T;xx=x;yy=y;zz=z;vx=v(l) ;vy=v(2) ;vz=v(3) ; 
else  xx= [xx  x] ;yy= [yy  y] ; zz= [zz  z] ; t= [t  ts+T] ; . . . 
vx=[vx  v ( 1 ) ] ; vy= [vy  v(2) ] ;vz=[vz  v(3) ] ; 


end; 


%Estimated  RDOA  error 
errt=norm(d  -  db) ; 

%Fill  matrix  for  estimated  RDOA  error 
if  ts==0;xet=errt; 
else  xet= [xet  errt] ; 
end; 

% Apply  Smith -Abel  algorithm  to  find  estimated  position  from  RDOAs 
source_est=smith_abel (RL,  d, w, R) ; 
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%Estimated  position  error 

errp=norm( sourcel  -  source_est) ; 

%Fill  matrix  for  estimated  position  error 
if  ts==0 ;xep=errp; 
else  xep=[xep  errp] ; 
end 

%Return  to  trajectory  loop 
end; 

%Obtain  position  and  velocity  components  vectors 
pos=[t'  xx'  yy'  zz'  ]  ; 
pos=pos ' ; 

v=  [  t '  vx '  vy '  vz ' ] ; 

V— v '  ; 

%Output  results  to  screen 
merrt=mean(xet) ; 

fprintf ( 'Mean  RDOA  Error  (ft)  =  %g  \n'  ,merrt) ; 
merrp=mean (xep) ; 

fprintf ( 'Mean  Position  Error  (ft)  =  %g  \n'  ,merrp) ; 

%Plot  graphs 

figure (1) 

plot (t,xet) ; grid; 

title ('RDOA  Error  vs  Time'); 

xlabel ( 'Time  -  sec ') ;ylabel ( 'RDOA  error  -  ft'); 

figure (2 ) 

plot (t,xep) ; grid; 

title (' Position  Error  vs  Time'); 

xlabel ('Time  -  sec ') ;ylabel (' Position  Error  -  ft'); 
figure (3) 

plot 3  (xx, yy,  zz)  ;grid; 
hold  on 

plot3 (RL( : ,1) ,RL( : ,2) ,RL( : ,3) , '*') ; 
hold  off 

title ( 'Missile  Trajectory  and  Remote  Receiver  Locations'); 
xlabel ( ' ft ' ) ;ylabel ( ' ft ' ) ; zlabel ( ' ft ' ) ; 


trajallmnc  .m 

%Position  from  FDOA  using  Smith^Abel  method 

%Inputs  missile  trajectory  and  receiver  SNR  and  calculates 

%RDOA  from  FDOA  and  position  with  SNR  measurement  error 

%Random  white  Gaussian  noise  added  to  signal 

%Monte  Carlo  simulation,  15-40  dB,  600  iterations 

%A11  receivers  utilized  throughout  entire  trajectory 

%trajallnmc  .m 

%Load  White  Sands  missile  trajectory 
load  traj ; 
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t=traj (1, : ) ;xm=traj (2, : ) ;ym=traj (3,  : ) ; zm=tra j (4,  : )  ; 

%Smooth  missile  trajectory 
xms=polyval (polyf it (t,xm, 14) ,  t)  ; 
yms=polyval (polyf it ( t ,ym, 14) ,  t) ; 
zms=polyval (polyf it ( t , zm, 14)  ,  t)  ; 

%Obtain  polynomial  coefficients  for  trajectory 
xpc=polyfit { t,xm, 14) ; 
ypc=polyf it ( t ,ym, 14) ; 
zpc=polyf it ( t , zm, 14) ; 

%Def ine  coefficients  for  velocity  (first  derivative  of  polynomials) 
i=l : 14 ; 

xvc= ( 1 5  —  i )  . *xpc (i) ; 
yvc=(15-i)  .  *ypc(i)-; 
zvc= (15-i) . *zpc (i) ; 

% Evaluate  velocity  polynomials  to  obtain  velocity  component  vectors 
vxs=polyval  (xvc,  t)  ; 
vys=polyval  (yvc,  t)  ; 
vzs=polyval  (zvc,  t)  ; 

%White  Sands  remote  receiver  positions  (feet) 

RL=[-9243.4  -3085.7  -127;  -9193.6  -3072.3  -125;  -9278.7  -3201.8 

-92.1;  -8992.9  -2879.8  -109.3;  1757  430.4  40.4;  -6458.5  -2112.3 
-39.3;  -3704.2  -1214.7  -33.6;  607  -3517.8  208.3;  -5305.4  2653 
-153.6] ; 

%Initialize  range  difference  vectors 

d=zeros (8,1) ; 

da=zeros (8,1) ; 

di=zeros (8,1) ; 

db=zeros (8,1) ; 

%Range  from  non-reference  receiver  to  reference  receiver 
R= zeros (8,1) ; 

R ( 1 )  =  norm(RL ( 1 , : )  -  RL(2,:)); 

R (2 )  =  norm ( RL ( 1 , : )  -  RL(3,:)); 

R (3 )  =  norm(RL(l, :)  -RL(4,:)); 

R ( 4 )  =  norm ( RL ( 1 , : )  -  RL ( 5 , : ) ) ; 

R ( 5 )  =  norm(RL(l,:)  -  RL(6,:)); 

R ( 6 )  =  norm ( RL ( 1 , : )  -  RL(7,:)); 

R ( 7 )  =  norm (RL ( 1 , : )  -  RL(8,:)); 

R ( 8 )  =  norm(RL(l,:)  -  RL(9,:)); 

%Set  all  weights  equal  to  one  (using  all  receivers) 
w=ones (8,1) ; 

%Telemetry  and  sampling  data 


f o=2 . 3e9 ;  %  Hz 

wl=(3e8/fo) /.3048;  %  feet 

T=le-2 ;  %  sec 

Tmax=6.41-;  %  sec 


%Start  S3STR  loop  for  receiver  SNR  values  (15-40  dB) 
for  snrdb=15:40; 
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snr=10A (snrdb/10)  ; 


%Theoretical  minimum  frequency  RMS  measurement  error 
dfrms=sqrt (3/snr) / (pi*T) ; 

%Start  Monte  Carlo  loop  (600  iterations) 
for  count=l:600; 


%Initial  range  differences,  with  Receiver  1  as  the  reference 
sourcei= [xms ( 1)  yms(l)  zms(l)]; 


di(l) 
di  (2 ) 
did) 
di  (4) 
did) 
di  (6) 
di  (7 ) 
di  (8) 


norm) sourcei 
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))  - 
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))  - 
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)  ) 
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d=di ; 


%Generate  loop  for  sampling  trajectory  points 

%Read  positions  and  velocities  at  each  sample  point  (every  .01  sec) 
for  t s  =  0 : T : Tmax ; 
sample=ts/T+l; 

sample=round ( sample)  ;  %Sample  must  be  an  integer  value 

x=xms  (1,  sample)  ;y=yms  (1,  sample)  ;  z=zms ( 1 , sample)  ; 

xl=xms  (1,  sample+1)  ;yl=yms  (1 ,  sample+1 )  ;  zl-zms  (1 ,  sample+1)  ; 

v=[vxs(l, sample)  vys  ( 1 ,  sample)  vzs (1, sample) ] ; 

vl= tvxs (1, sample+1)  vys (1, sample+1)  vzs (1 , sample+1)]; 


%Find  source  position  at  present  time 
source= [x  y  z] ; 


%Find  source  position  at  next  time 
sourcel=[xl  yl  zl]  ; 

%Find  source-to-receiver  vectors  at  present  time 
Bl=source-RL (1,:); 

B2=source-RL(2, : ) ; 

B3=source-RL (3 ,  : )  ; 

B4=source-RL (4 , : ) ; 

B5-source-RL ( 5 , : ) ; 

B6=source-RL ( 6 , : ) ; 

B7 =source-RL (7,:); 

B8=source-RL(8/ : ) ; 

B9=source-RL(9, : ) ; 


%Find  Doppler  frequency  vector  at  present  time 
F=v/wl ; 

%Find  Doppler  frequencies  at  receivers  at  present  time  (with  noise) 
fdl=dot  (Bl,  F)  /norm(Bl)  +df rms*randn; 
fd2=dot  (B2,F)  /norm(B2)  +df rms*randn; 
fd3=dot  (B3  ,  F)  /norm (B3  )  +dfrms*randn; 
fd4=dot  (B4 ,  F)  /norm(B4) +dfrms*randn; 
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fd5=dot (B5 , F) /norm(B5) +dfrms*randn; 
fd6=dot  (B6, F)  /norm(B6) +dfrms*randn; 
fd7=dot (B7 , F) /norm (B7 ) +dfrms*randn; 
fd8=dot (B8 ,  F) /norm(B8) +dfrms*randn; 
fd9=dot (B9 , F) /norm(B9) +df rms*randn; 

%Find  differential  RDOAs  at  present  time 
drdoa2a=-wl*T* ( fdl-f d2 ) ; 
drdoa3a=-wl*T* ( fdl-f d3 ) ; 
drdoa4a=-wl*T* (fdl-fd4) ; 
drdoa5a=-wl*T* ( fdl-f d5) ; 
drdoa6a=-wl*T* (fdl-fd6) ; 
drdoa7a=-wl*T* (fdl-fd7) ; 
drdoa8a=-wl*T* (fdl-fd8) ; 
drdoa9a=-wl*T* ( fdl-f d9 ) ; 

%Find  source-to-receiver  vectors  at  next  time 
Bl=sourcel-RL ( 1 # : ) ; 

B2=sourcel-RL (2 , : ) ; 

B3=sourcel-RL (3 , : ) ; 

B4=sourcel-RL (4, : ) ; 

B5=sourcel-RL ( 5 , : ) ; 

B6=sourcel-RL ( 6, : ) ; 

B7=sourcel-RL(7, : ) ; 

B8=sourcel-RL ( 8 , : ) ; 

B9=sourcel-RL (9, : ) ; 

%Find  Doppler  frequency  vector  at  next  time  ' 

Fl=vl/wl; 

%Find  Doppler  frequencies  at  receivers  at  next  time  (with  noise) 
fdl=dot (Bl#  FI) /norm(Bl) +dfrms*randn; 
fd2=dot (B2/F1) /norm  (B2 ) +df rms*randn; 
fd3=dot (B3 , FI) /norm(B3 ) +df rms*randn; 
fd4=dot (B4, FI) /norm(B4 ) +df rms*randn; 
fd5=dot (B5, FI) /norm(B5) +dfrms*randn; 
fd6=dot (B6, FI) /norm (B6 ) +df  rms*randn; 
fd7=dot (B7,F1) /norm  ( B7 ) +df rms  *randn ; 
fd8=dot (B8, FI) /norm(B8 ) +df rms*randn; 
fd9=dot (B9 ,  FI ) /norm{B9) +df rms*randn; 

%Find  differential  RDOAs  at  next  time 
drdoa2b=-wl*T* (fdl-f d2) ; 
drdoa3b=-wl*T* ( fdl-f d3) ; 
drdoa4b=-wl*T* (fdl-fd4) ; 
drdoa5b=-wl*T* (fdl-fd5) ; 
drdoa6b=-wl*T* (fdl-fd6) ; 
drdoa7b=-wl*T* (fdl-fd7) ; 
drdoa8b=-wl*T* (fdl-fd8) ; 
drdoa9b=-wl*T* ( fdl-fd9 ) ; 

%Find  average  differential  RDOAs 

drdoa2= (drdoa2a+drdoa2b) /2 ; 
drdoa3= (drdoa3a+drdoa3b) /2 ; 
drdoa4= (drdoa4a+drdoa4b) /2 ; 
drdoa5= (drdoa5a+drdoa5b) /2 ; 
drdoa6= (drdoa6a+drdoa6b) / 2 ; 
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drdoa7= (drdoa7a+drdoa7b) / 2 ; 
drdoa8= (drdoa8a+drdoa8b) / 2 ; 
drdoa9= (drdoa9a+drdoa9b) /2  ; 

%Find  RDOAs 

da ( 1 ) =d ( 1 ) +drdoa2  ; 
da ( 2 ) =d ( 2 ) +drdoa3  ; 
da  ( 3 ) =d ( 3 ) +drdoa4  ; 
da ( 4 ) =d ( 4 ) +drdoa5  ; 
da ( 5 ) =d ( 5 ) +drdoa6  ; 
da ( 6 ) =d ( 6 ) +drdoa7  ; 
da  ( 7 ) =d ( 7 ) +drdoa8  ; 
da ( 8 ) =d ( 8 ) +drdoa9  ; 

%Reset  range  differences 
d=da; 

%Find  RDOAs  from  geometry 

db(l)  =  norm(sourcel  -  RL(2,:))  -  norm(sourcel  -  RL(1,:)); 

db (2 )  =  norm(sourcel  -  RL(3,:))  -  norm(sourcel  -  RL(1,:)); 

db (3 )  =  norm(sourcel  -  RL(4,:))  -  norm(sourcel  -  RL(1,:)); 

db(4)  =  norm(sourcel  -  RL(5,:))  -  norm(sourcel  -  RL(1,:)); 

db(5)  =  norm(sourcel  -  RL(6,:))  -  norm(sourcel  -  RL(1,:)); 

db(6)  =  norm(sourcel  -  RL(7,:))  -  nonn(sourcel  -  RL ( 1 , : ) ) ; 

db (7)  =  norm(sourcel  -  RL(8,:))  -  norm(sourcel  -  RL(1,:)); 

db(8)  =  norm(sourcel  -  RL(9,:))  -  norm(sourcel  -  RL(1,:)); 

%Estimated  RDOA  error 

errt=norm(d  -  db) ; 

%Fill  matrix  for  estimated  RDOA  error 
if  ts==0 ;xet=errt; 
else  xet= [xet  errt] ; 
end; 

% Apply  Smith -Abel  algorithm  to  find  estimated  position  from  RDOAs 
source_est=smith_abel (RL,d,w,R) ; 

%Estimated  position  error 

errp=norm (sourcel  -  source_est) ; 

%Fill  matrix  for  estimated  position  error 
if  ts==0 ;xep=errp; 
else  xep= [xep  errp] ; 
end 

%Return  to  trajectory  loop 
end; 

%Calculate  mean  RDOA  and  position  errors  for  missile  trajectory 

%Fill  matrices  for  these  values 
merrt=mean (xet) ; 
merrp=mean (xep) ; 

if  count==l ;meanxet=merrt;meanxep=merrp; 

else  meanxet= [meanxet  merrt] ;meanxep= [meanxep  merrp] ;end 


%Return  to  Monte  Carlo  loop 
end; 

%Calculate  and  print  mean  values  for  Monte  Carlo  iterations 
mcet=mean (meanxet ) ; 
mce=mean (meanxep) ; 
f print f ( ' SNR (dB)  =  %g\n',  snrdb); 
fprintf (' Error  in  RDOA  (ft)  =  %g\n',  meet); 
fprintf ( 'Error  in  Position  (ft)  =  %g\n'/  mce) ; 

%Fill  SNR,  RDOA,  and  position  error  matrices  for  plotting 
if  snrdb==15 ;  snrdbm=  snrdb; me etm=mcet ;mcem=mce; 

else  snrdbm=  [ snrdbm  snrdb]  ; mce tm=  [mcetm  meet ]  ;mcem=  [mcem  mce] ; end; 

%Return  to  SNR  loop 
end; 

%Save  data  to  workspace 
save  data  snrdbm  mcetm  mcem; 

%Plot  graphs 
figure (1) 

plot (snrdbm, mcetm) ;grid; 

title ('Range  Difference  Error  vs.  SNR'); 

xlabel('SNR  -  dB ') ;ylabel (' Range  Difference  Error  -  ft'); 
figure (2 ) 

plot (snrdbm, mcem) ;grid; 

title (' Position  Error  vs.  SNR'); 

xlabel ( ' SNR  -  dB' ) ;ylabel (' Position  Error  -  ft'); 


angles,  m 

%Calculates  azimuth  and  elevation  angles  for  each  receiver  vs.  time 
%The  matrix  "pos"  obtained  from  executing  trajall.m  or  trajalln.m 
%angles .m 

%Load  trajectory  data  (4x642  matrix) 

%First  row  is  time,  other  rows  are  x,y,z  positions 
load  pos; 

%Define  receiver  location  vectors 
bsla= [ -9243 . 4  -3085.7  -127.0]; 
bslb=[-9193. 6  -3072.3  -125.0]; 
bslp= [ -9278 . 7  -3201.8  -92.1]; 
bsl0= [ -8992 . 9  -2879.8  -109.3]; 
bsll= [1757.0  430.4  40.4] ; 
bsl2= [-6458 . 5  -2112.3  -39.3]; 
bsl4= [ -3704 . 2  -1214.7  -33.6]; 
bsl5= [607 . 0  -3517.8  208.3]; 
bsl6= [ -5305 . 4  2653.0  -153.6]; 

%Translate  trajectory  such  that  receivers  are  at  origin 
trajla= [pos ( 2 , : ) -bsla (1) ;pos (3 , : ) -bsla (2) ;pos (4 , : ) -bsla (3 ) ] ; 
trajlb= [pos (2 , : ) -bslb (1) ;pos (3 , : ) -bslb (2) ;pos (4 , : ) -bslb (3 ) ] ; 
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trajlp= [pos (2 , :)-bslp(l) ;pos(3, :)-bslp(2) ;pos(4, :)-bslp(3)] ; 
t raj  10= [pos ( 2  ,  :  ) -bsl 0 ( 1 ) ; pos (3 , : ) -bslO { 2 ) ; pos ( 4 ,  : ) -bslO ( 3 ) ] ; 
traj 11= [pos (2 , : ) -bsll(l) ;pos(3, : ) -bsll(2) ;pos(4, : ) -bsll (3) ] ; 
traj 12= [pos (2 , :)-bs!2(l) ;pos(3, :)-bsl2(2) ;pos(4, :)-bsl2(3)] ; 
traj 14= [pos (2 ,  :  ) -bsl4 (1) ;pos (3 ,  : ) -bsl 4 (2 ) ;pos (4, : ) -bsl 4  (3 ) ] ; 
traj 15= [pos (2 , :)-bsl5(l) ;pos(3, :)-bsl5(2) ;pos(4, :)-bsl5(3)] ; 
traj 16= [pos ( 2 , : ) -bsl 6 (1) ;pos ( 3 , : ) -bsl 6 ( 2 ) ; pos ( 4 , : ) -bsl 6 ( 3 ) ] ; 

%Calculate  azimuth  (theta)  and  elevation  (phi)  angles 
%Convert  radians  to  degrees 

[ thetala  phila  rla] =cart2sph(trajla (1, :) ,trajla(2, :) /trajla(3/ :) ) ; 
thetala=thetala*180/pi ; 
phila=phila*180/pi ; 

[thetalb  philb  rib] =cart2sph(trajlb(l/ : ) , traj lb (2, : ) , traj lb (3, : ) ) ; 
thetalb=thetalb*180/pi ; 
philb=philb*180/pi ; 

[thetalp  philp  rip] =cart2sph (trajlp (1, : ) , trajlp (2 , : ) , trajlp  (3 ,  : ) )  ; 

thetalp=thetalp*180/pi; 

philp=philp*180/pi ; 

[ thetalO  philO  rlO] =cart2sph ( trajlO (1, : ) , traj 10 (2, : ) , traj 10 (3,:)); 
thetal0=thetal 0*180 /pi ; 
phi 10=phil 0*180 /pi ; 

[ thetall  phi  11  rll]  =cart2sph ( traj  11  (1,  :)  /trajll(2/  :)  .trajllO,  :)); 

thetall=thetall*180/pi ? 

phill=phill*180/pi; 

[thetal2  phil2  rl2 ] =cart2sph ( traj 12 (1, : ) , traj 12 (2, : ) , traj 12 (3, : ) ) ; 
thetal2=theta!2*180/pi ; 
phil2=phil2  *180/pi ; 

[thetal4  phil4  rl4] =cart2sph(trajl4 (1, :) /trajl4(2/ :) /trajl4(3/ :) ) ; 
thetal4=thetal4*180/pi ; 
phil4=phil4*180/pi ; 

[thetalS  phi 15  rl5] =cart2sph ( trajlS (1, :) /trajl5(2/ :) ,trajl5(3, :) ) ; 

thetal5=thetal5*180/pi ; 

phil5=phil5*180/pi; 

[thetal6  phil6  rl6] =cart2sph ( traj 16 (1, :) /trajl6(2/ :) ,trajl6(3, :) ) ; 
thetal6=thetal6*180/pi ? 
phil6=phil6*180/pi ; 

%Plot  graphs  of  theta  and  phi  versus  trajectory  time 
figure (1) ; 

plot (pos (1, : ) , thetala) ; 
hold  on; 

plot (pos (1, : ) / phila, ' r' ) ;grid; 
hold  off; 

title ( 'Angles  for  Receiver  la  -  Reference  Receiver'); 
xlabeM'Time  -  sec ');  ylabel  ( 'Angle  -  deg'); 

figure (2 ) ; 

plot (pos ( 1 , : ) , thetalb) ; 
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hold  on; 

plot (pos (1, : ) , philb, 'r ' ) ;grid; 
hold  off; 

title ( 'Angles  for  Receiver  lb'); 

xlabel ( 'Time  -  sec ') ;ylabel ( 'Angle  -  deg') 

figure (3) ; 

plot (pos (1, : ) , thetalp) ; 
hold  on; 

plot (pos (1, : ) ,philp, 'r' ) ;grid; 
hold  off; 

title ( 'Angles  for  Receiver  Ip'); 

xlabel ( 'Time  -  sec ') ;ylabel ( 'Angle  -  deg') 

figure (4) ; 

plot (pos (1, : ) , thetalO) ; 
hold  on; 

plot (pos (1/ : ) , philO, 'r ' ) ;grid; 
hold  off; 

title ( 'Angles  for  Receiver  10'); 

xlabel ('Time  -  sec ') ;ylabel ( 'Angle  -  deg') 

figure (5)  ; 

plot (pos (1, : )  , thetall) ; 
hold  on; 

plot (pos (1, : ) ,phill, 'r' ) ;grid; 
hold  off; 

title { 'Angles  for  Receiver  11'); 

xlabel {'Time  -  sec ') ;ylabel ( 'Angle  -  deg') 

figure (6) ; 

plot (pos ( 1 , : ) / thetal2 ) ; 
hold  on; 

plot (pos ( 1 , : ) , phil2 , ' r ' ) ; grid; 
hold  off; 

title ( 'Angles  for  Receiver  12'); 

xlabel ('Time  -  sec ') ;ylabel ( 'Angle  -  deg') 

figure { 7 ) ; 

plot (pos  ( 1 ,  : ) y  thetal4 ) ; 
hold  on; 

plot (pos (1, : ) ,phil4, 'r' ) ;grid; 
hold  off; 

title ( 'Angles  for  Receiver  14'); 

xlabel ('Time  -  sec ') ;ylabel ( 'Angle  -  deg') 

figure { 8 ) ; 

plot (pos ( 1 , : ) , thetalS ) ; 
hold  on; 

plot (pos (1, : ) ,phil5, 'r' ) ;grid; 
hold  off; 

title ( 'Angles  for  Receiver  15'); 

xlabel ('Time  -  sec ') ;ylabel ( 'Angle  -  deg') 

figure (9) ; 

plot (pos (1, : ) , theta!6) ; 
hold  on; 
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plot (pos (1, : ) ,phil6, 'r ' ) ;grid; 
hold  off; 

title ( 'Angles  for  Receiver  16'); 

xlabel ( ' Time  -  sec ') ;ylabel ( 'Angle  -  deg'); 


t  raj  group,  m 

%Position  from  FDOA  using  Smith-Abel  algorithm 
%Inputs  missile  trajectory  and  receiver  SNR  and  calculates 
%RDOA  from  FDOA  and  position  with  SNR  measurement  error 
%Only  receivers  with  telemetry  in  3-dB  beam  width  utilized 
%trajgroup.m 

%Load  White  Sands  missile  trajectory 
load  tra j ; 

t=traj {1, : ) ;xm=traj (2, : ) ;ym=traj (3, : ) ; zm=traj (4,  :  )  ; 

%Smooth  missile  trajectory 
xms=polyval (polyf it { t ,  xm, 14)  ,  t)  ,* 
yms=polyval {polyf it ( t ,ym, 14) , t ) ; 
zms=polyval (polyf it ( t , zm, 14 ) , t ) ; 

%Obtain  polynomial  coefficients  for  trajectory 
xpc=polyf it (t,xm, 14) ; 
ypc=polyf it { t ,ym, 14) ; 
zpc=polyfit (t, zm, 14) ; 

%Def ine  coefficients  for  velocity  (first  derivative  of  polynomials) 
i  =  1 : 1 4 ; 

xvc= (15-i) . *xpc (i) ; 
yvc= (15-i) . *ypc (i) ; 
zvc= (15-i)  . *zpc  (i) ; 

%Evaluate  velocity  polynomials  to  obtain  velocity  component  vectors 
vxs=polyval (xvc, t) ; 
vys=polyval (yvc , t ) ; 
vzs=polyval (zvc, t) ; 

%White  Sands  remote  receiver  positions  (feet) 

RL= [-9243.4  -3085.7  -127;  -9193.6  -3072.3  -125;  -9278.7  -3201.8 

-92.1;  -8992.9  -2879.8  -109.3;  1757  430.4  40.4;  -6458.5  -2112.3 
-39.3;  -3704.2  -1214.7  -33.6;  607  -3517.8  208.3;  -5305.4  2653 
-153.6]  ; 

%Initialize  range  difference  vectors 

d=zeros (8,1)  ; 

da=zeros (8,1) ; 

di=zeros (8,1); 

db=zeros (8,1) ; 

%Range  from  non-reference  receiver  to  reference  receiver 
R=zeros (8,1); 

R ( 1 )  =  norm ( RL ( 1 , : )  -  RL ( 2 , : ) ) ; 

R ( 2 )  =  norm(RL(l, : )  -RL(3,:)); 

R(3)  =  norm(RL ( 1 , : )  -  RL(4,:)); 
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R ( 4 )  =  norm ( RL ( 1 ,  : )  -  RL ( 5 , : ) ) ; 

R ( 5 )  =  norm ( RL ( 1 , : )  -  RL ( 6 , : ) ) ; 

R ( 6 )  =  norm ( RL ( 1 , : )  -  RL ( 7 , : ) ) ; 

R (7)  =  norm (RL ( 1 , : )  -  RL(8,:)); 

R ( 8 )  =  norm(RL(l,:)  -  RL(9,:)); 

%Telemetry  and  sampling  data 
f o=2 . 3e9 ;  %  Hz 

wl=(3e8/fo) /.3048;  %  feet 

T=le-2 ;  %  sec 

Tmax=6.41;  %  sec 

%Initial  range  differences,  with  Receiver  1  as  the  reference 
sourcei=  [xms  (1)  yms(l)  zms(l)]; 

di(l)  =  norm(sourcei  -  RL(2,:))  -  norm(sourcei  -  RLfl,:)); 

di ( 2 )  =  norm(sourcei  -  RL(3,:))  -  norm(sourcei  -  RL(1,:)); 

di { 3 )  =  norm(sourcei  -  RL{4,:))  -  norm ( sourcei  -  RLfl,:)); 

di(4)  =  normfsourcei  —  RL ( 5 , : ) )  -  norm(sourcei  -  RLfl,:)); 

di (5 )  =  norm(sourcei  -  RL ( 6 , : ) )  -  norm(sourcei  -  RL(1,:)); 

di ( 6 )  =  norm(sourcei  -  RL(7,:))  -  normfsourcei  -  RLfl,:)); 

di (7 )  =  norm(sourcei  -  RL(8,:))  -  normfsourcei  -  RL ( 1 , : ) ) ; 

di ( 8 )  =  normfsourcei  -  RL ( 9 # : ) )  -  normfsourcei  -  RL ( 1 , : ) ) ; 

d=di  ; 

% Input  receiver  SNR 
snrdb= input { ' snr (dB)  = ' )  ; 
snr=10A ( snrdb/10 ) ; 

%Theoretical  minimum  frequency  RMS  measurement  error 
dfrms=sqrt (3/snr) / (pi*T) ; 

%Generate  loop  for  sampling  trajectory  points 

%Read  positions  and  velocities  at  each  sample  point  (every  .01  sec) 
for  ts=0:T:Tmax; 
sample=ts/T+l; 

sample=round ( sample)  ;  %Sample  must  be  an  integer  value 

x=xms  (1,  sample)  ;y=yms  (1,  sample)  ;  z=zms  ( 1 ,  sample)  ; 

xl=xms  (1,  sample+1)  ;yl=yms  (1,  sample+1)  ;  zl=zms (1, sample+1)  ; 

v=  [ vxs  (1,  sample)  vys  (1 ,  sample)  vzs  fl,  sample)  ]  ; 

vl= [vxs (1, sample+1)  vys ( 1 , sample+1 )  vzs ( 1 , sample+1) ] ; 

%Find  source  position  at  present  time 
source= [x  y  z] ; 

%Find  source  position  at  next  time 
sourcel=[xl  yl  zl] ; 

%Find  source-to-receiver  vectors  at  present  time 
Bl=source-RL (1,  : )  ; 

B2=source-RL (2 ,  : )  ; 

B3=source-RL (3 , : ) ; 

B4=source-RL ( 4 ,  : )  ; 

B5=source-RL ( 5 , : ) ; 

B6=source-RL(6,:); 

B7=source-’RL(7,  :  )  ; 
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B8=source-RL ( 8 , : ) ; 

B9=source-RL ( 9  ,  : )  ; 

%Find  Doppler  frequency  vector  at  present  time 
F=v/wl; 

%Find  Doppler  frequencies  at  receivers  at  present  time 
f dl=dot (B1 , F) /norm { B1 ) +df rms; 
f d2=dot (B2 ,  F) /norm(B2) +dfrms; 
fd3=dot (B3 , F) /norm(B3 ) +dfrms; 
fd4=dot (B4, F) /norm ( B4 ) +df  rms ; 
fd5=dot (B5 / F) /norm(B5) +dfrms; 
fd6=dot (B6, F) /norm{B6) +dfrms; 
fd7=dot (B7 ,  F) /norm ( B7 ) +df rms; 
f d8=dot (B8 , F) /norm {B8 ) +df rms ; 
fd9=dot (B9, F) /norm(B9) +dfrms; 

%Find  differential  RDOAs  at  present  time 
drdoa2a=-wl*T* (fdl-fd2) ; 
drdoa3a=-wl*T* (fdl-fd3) ; 
drdoa4a=-wl*T* (fdl-fd4) ; 
drdoa5a=-wl*T* (fdl-fd5) ; 
drdoa6a=-wl*T* (fdl-fd6)  ; 
drdoa7a=-wl*T* (fdl~fd7) ; 
drdoa8a=-wl*T* (fdl-fd8) ; 
drdoa9a=-wl*T* (fdl-fd9) ; 

%Find  source-to-receiver  vectors  at  next  time 
Bl  =  sourcel~RL(l/  : )  ; 

B2=sourcel-RL (2, : ) ; 

B3=sourcel-RL (3  ,  : )  ; 

B4=sourcel-RL (4, : ) ; 

B5  =  sourcel-RL (5,  : )  ; 

B6=sourcel-RL(6,  : )  ; 

B7  =  sourcel-RL(7/  :  )  ; 

B8  =  sourcel-RL ( 8 ,  :  )  ; 

B9=sourcel-RL(9/ : ) ; 

%Find  Doppler  frequency  vector  at  next  time 
Fl=vl/wl; 

%Find  Doppler  frequencies  at  receivers  at  next  time 
fdl=dot (B1 , FI ) /norm ( B1 ) +df rms ; 
fd2=dot (B2,F1) /norm (B2 ) +dfrms; 
fd3=dot (B3 , FI) /norm(B3) +dfrms; 
f d4=dot (B4 , FI) /norm ( B4 ) +df rms ; 
f d5=dot ( B5 , FI ) /norm ( B5 ) +df rms ; 
fd6=dot (B6,F1) /norm(B6) +dfrms; 
fd7=dot (B7 , FI) /norm(B7) +dfrms ; 
fd8=dot (B8 , FI ) /norm (B8) +dfrms; 
f d9=dot (B9 , FI ) /norm  (B9 ) +df rms; 

%Find  differential  RDOAs  at  next  time 
drdoa2b=-wl*T* (fdl-fd2) ;  v 

drdoa3b=-wl*T* (fdl-fd3) ; 
drdoa4b=-wl*T* (fdl-fd4) ; 
drdoa5b=-wl*T* (fdl-fd5) ; 
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drdoa6b=-wl*T* (fdl-fd6) ; 
drdoa7b=-wl*T* ( f dl-fd7 ) ; 
drdoa8b=-wl*T* ( f dl-f d8 ) ; 
drdoa9b=-wl*T* (fdl-fd9) ; 

%Find  average  differential  RDOAs 
drdoa2= (drdoa2a+drdoa2b) /2 ; 
drdoa3= (drdoa3a+drdoa3b) /2; 
drdoa4= (drdoa4a+drdoa4b) / 2; 
drdoa5= (drdoa5a+drdoa5b) /2 ; 
drdoa6= (drdoa6a+drdoa6b) /2; 
drdoa7= (drdoa7a+drdoa7b) /2; 
drdoa8= (drdoa8a+drdoa8b) /2; 
drdoa9= (drdoa9a+drdoa9b) /2; 

%Find  RDOAs 

da ( 1 ) =d ( 1 ) +drdoa2  ; 
da ( 2 ) =d ( 2 ) +drdoa3 ; 
da ( 3 ) =d ( 3 ) +drdoa4 ; 
da ( 4 ) =d ( 4 ) +drdoa5 ; 
da ( 5 ) =d ( 5 ) +drdoa6 ; 
da ( 6 ) =d ( 6 ) +drdoa7 ; 
da ( 7 ) =d ( 7 ) +drdoa8 ; 
da ( 8 ) =d ( 8 ) +drdoa9 ; 

%Reset  range  differences 
d=da  ; 

%Find  RDOAs  from  geometry 

db(l)  =  normtsourcel  -  RL(2,:))  -  normtsourcel  -  RL(1,:)); 

db(2)  =  norm(sourcel  -  RL(3,:))  -  norm(sourcel  -  RL(1,:)); 

db(3)  =  norm(sourcel  -  RL(4,:))  -  norm(sourcel  -  RL(1,:)); 

db(4)  =  normtsourcel  -  RL(5,:))  -  normtsourcel  -  RL(1,:)); 

db(5)  =  normtsourcel  -  RL(6,:))  -  normtsourcel  -  RL(1,:)); 

db(6)  =  normtsourcel  -  RL(7,:))  -  normtsourcel  -  RL(1,:)); 

db(7)  =  normtsourcel  -  RL(8,:))  -  normtsourcel  -  RL(1,:)); 

db(8)  =  normtsourcel  -  RL(9,:))  -  normtsourcel  -  RL(1,:)); 

%Fill  matrices  for  time,  position,  and  velocity  components 
if  ts==0 ; t=ts+T; xx=x,-yy=y ;  zz=z ; vx=v ( 1 )  ;vy=v(2)  ;vz=v(3)  ; 
else  xx= [xx  x] ;yy= [yy  y] ; zz= [zz  z];t=[t  ts+T] ; . . . 
vx— [vx  v(l) ] ; vy= [vy  v(2) ] ;vz= [vz  v(3) ] ; 

end; 

%Use  Smith-Abel  algorithm  to  obtain  estimated  position  from  RDOAs 

%Apply  45  deg  3-dB  beam  width  of  receiver  antennas  for  signal  reception 

%Weights  of  1  define  receivers  used,  weights  of  0  define  those  not  used 

%Receivers  la, lb, lp, 10, 11, 12,14  and  15  within  view  at  time  start 
if  ts<2.12 

w= [ 1 ; 1 ; 1 ; 1 ; 1 ; 1 ; 1 ; 0 ] ; 

du=[d(l) ; d ( 2 ) ; d ( 3 ) ; d ( 4 ) ; d ( 5 ) ; d ( 6 ) ;  d  ( 7 ) ] ; 
dbn=[db(l) ;db{2) ;db(3) ;db(4) ;db(5) ;db(6) ;db(7) ] ; 
source_est=smith_abel (RL,d,w,R) ; 
end; 

%Lose  Receiver  12  at  2 . 12  sec 
if  ts>=2 . 12  Sc  ts<3 . 29 
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w= [ 1 ; 1 ; 1 ; 1 ; 0 ; 1 ; 1 ; 0 ]  ; 
du= [ d { 1 ) ; d ( 2 ) ; d { 3 ) ; d ( 4 ) ; d ( 6 ) ; d ( 7 ) ] ; 
dbn= [ db ( 1 ) ;db(2) ;db(3) ;db(4) ;db(6) ;db(7) ] ; 
source_est=smith_abel ( RL ,  d ,  w , R) ; 
end; 

%Lose  Receiver  14  at  3.29  sec 
if  ts>=3.29  &  ts<4.49 

W=  [  1 ; 1 ; 1 ; 1 ; 0 ; 0 ; 1 ;  0  ]  ; 
du= [ d ( 1 ) ; d ( 2 ) ; d ( 3 ) ; d ( 4 ) ; d ( 7 ) ] ; 
dbn=  t  db ( 1 ) ;db(2);db(3);db(4);db(7)]; 
source_est=smith_abel (RL,  d,  w,  R) ; 
end; 

%Obtain  RDOA  for  added  receiver  from  estimated  position 
if  ts==4 . 49 

d ( 8 ) =norm ( source_es t  -  RL(9,:))  -  norm( source^est  -  RL ( 1 , : ) )  ; 
d(8) =d(8) +drdoa9; 
end; 

%Lose  Receiver  15  and  gain  Receiver  16  at  4.49  sec 
if  ts>=4.49  &  ts<5.95 
w=  [  1 ;  1 ;  1 ;  1 ;  0  ;  0 ;  0 ;  1  ]  ; 
du=  [  d  ( 1 )  ;  d  ( 2 )  ;  d  ( 3 )  ;  d  ( 4 )  ;  d  ( 8 )  ]  ; 
dbn= [db (1) ;db(2);db(3);db(4);db(8)]; 
source_est=smith_abel (RL, d, w, R) ; 
end; 

%Lose  Receiver  11  at  5.95  sec 
if  ts>=5 . 95  &  ts<=Tmax 
w=  [  1 ; 1 ; 1 ; 0 ; 0 ; 0 ; 0 ;  1  ]  ; 
du= [d ( 1) ; d ( 2 ) ;d(3) ;d(8)  ]  ; 
dbn=[db(l) ;db(2) ;db(3) ;db{8) ] ; 
source_est=smith_abel (RL/d^W/R) ; 
end; 

%Estimated  RDOA  error  of  used  range  differences 
errt=norm(du  -  dbn) ; 

%Fill  matrix  for  estimated  RDOA  error 
if  ts==0 ;xet=errt ; 
else  xet=[xet  errt] ; 
end; 

%Estimated  position  error 

errp=norm(sourcel  -  source_est) ; 

%Fill  matrix  for  estimated  position  error 
if  ts==0 ;xep=errp; 
else  xep=[xep  errp] ; 
end 

%Return  to  trajectory  loop 

end; 

%Obtain  position  and  velocity  components  vectors 

pos=[t'  xx'  yy'  zz' ] ; 
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pos=pos ' ; 

v= [ t 7  vx7  vy '  vz  7 ] ; 
v=v 7  ; 

%Output  results  to  screen 
merrt=mean (xet) ; 

fprintf ( 'Mean  RDOA  Error  (ft)  =  %g  \n7  ,merrt) ; 
merrp=mean (xep) ; 

fprintf ( 'Mean  Position  Error  (ft)  =  %g  \n7  ,merrp) ; 

%Plot  graphs 

figured) 

plot (t/Xet) ;grid; 

title ('RDOA  Error  vs  Time'); 

xlabel ( 'Time  -  sec7 ) ;ylabel ( 'RDOA  error  -  ft7); 

figure (2) 

plot (t, xep) ; grid; 

title ( 7 Position  Error  vs  Time7); 

xlabel ('Time  -  sec 7 ) ;ylabel { 7 Position  Error  -  ft'); 
figure (3 ) 

plot3 (xx, yy, zz ) ; grid; 
hold  on 

plot3 (RL( : ,1) ,RL( : ,2) ,  RL  ( : ,3) , ) ; 
hold  off 

title ( 'Missile  Trajectory  and  Remote  Receiver  Locations'); 
xlabel ( 7  f t 7 ) ;ylabel ( 7  ft 7 ) ; zlabel ( 7  f t 7 ) ; 


t ra j  groupn . m 

%Position  from  FDOA  using  Smith-Abel  algorithm 

%Inputs  missile  trajectory  and  receiver  SNR  and  calculates 

%RDOA  from  FDOA  and  position  with  SNR  measurement  error 

%Random  white.  Gaussian  noise  added  to  signal 

%Only  receivers  with  telemetry  in  3-dB  beam  width  utilized 

%  tra j  groupn . m 

%Load  White  Sands  missile  trajectory 
load  traj ; 

t=traj ( 1 ,  : ) ;xm=traj (2 , : )  ;ym-traj (3  7  : )  ; zm=traj (4, : ) ; 

%Smooth  missile  trajectory 
xms=polyval (polyfit (t.xm, 14) , t) ; 
yms=polyval (polyfit (t/ym, 14) , t) ; 
zms=polyval (polyfit (t , zm, 14) , t) ; 

%Obtain  polynomial  coefficients  for  trajectory 
xpc=polyf it ( t,xm/ 14); 
ypc=polyf it ( t ,ym/ 14)  ; 
zpc=polyf it (t , zm, 14)  ; 

%Define  coefficients  for  velocity  (first  derivative  of  polynomials) 
i=l : 14 ; 

xvc= (15-i) . *xpc (i) ; 
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yvc= ( 15  —  i ) . *ypc (i) ; 
zvc=(15-i) .*zpc(i) ; 

%Evaluate  velocity  polynomials  to  obtain  velocity  component  vectors 
vxs=polyval (xvc,  t)  ; 
vys=polyval (yvc, t) ; 
vzs=polyval (zvc, t)  ; 

%White  Sands  remote  receiver  positions  (feet) 

RL= [ -9243 . 4  -3085.7  -127;  -9193.6  -3072.3  -125;  -9278.7  -3201.8 

-92.1;  -8992.9  -2879.8  -109.3;  1757  430.4  40.4;  -6458.5  -2112.3 
-39.3;  -3704.2  -1214.7  -33.6;  607  -3517.8  208.3;  -5305.4  2653 
-153.6] ; 

%Initialize  range  difference  vectors 

d=zeros (8,1) ; 

da=zeros (8,1); 

di=zeros (8,1); 

db=zeros (8,1) ; 

%Range  from  non-reference  receiver  to  reference  receiver 
R=zeros (8,1) ; 

R ( 1 )  =  norm (RL ( 1 , : )  -  RL(2,:)); 

R ( 2 )  =  norm ( RL ( 1 , : )  -  RL ( 3 , : )  )  ; 

R  ( 3 )  =  norm ( RL ( 1 , : )  -  RL ( 4 ,  : ) ) ; 

R ( 4 )  =  norm ( RL ( 1 , : )  -  RL ( 5 , : ) ) ; 

R ( 5 )  =  norm(RL(l, :)  -RL(6,:)); 

R ( 6 )  =  norm ( RL ( 1 , : )  -  RL ( 7 , : ) ) ; 

R ( 7 )  =  norm ( RL ( 1 , : )  -  RL ( 8 , : ) ) ; 

R(8)  =  norm(RL(l, : )  -  RL ( 9 , : ) ) ; 

%Telemetry  and  sampling  data 
fo=2 . 3e9 ;  %  Hz 

wl=(3e8/fo) /.3048;  %feet 

T=le-2;  %  sec 

Tmax=6 . 41 ;  %  sec 

%Initial  range  differences,  with  Receiver  1  as  the  reference 
sourcei= [xms (1)  yms(l)  zms(l)]; 

did)  =  norm(sourcei  —  RL  (2 ,  ; ) )  -  norm(sourcei  -  RL(1,:)); 

di (2)  =  norm (sourcei  -  RL(3,:))  -  norm(sourcei  -  RL(1,:)); 

di  (3)  =  norm  (sourcei  -  RL(4,:))  -  norm(sourcei  -  RL(1,:)); 

di (4)  =  norm (sourcei  -  RL(5,:))  -  norm (sourcei  -  RL(1,:)); 

di ( 5 )  =  norm (sourcei  -  RL(6,:))  -  norm (sourcei  -  RL{1,:)); 

di ( 6 )  =  norm (sourcei  -  RL(7,:))  -  norm(sourcei  -  RL(1,:)); 

di (7)  =  norm (sourcei  -  RL(8,:))  -  norm (sourcei  -  RL(1,:)); 

di ( 8 )  =  norm (sourcei  -  RL(9,:))  -  norm (sourcei  -  RL(1,:)); 

d=di  ; 

%Input  receiver  SNR 
snrdb= input ( ' snr (dB) = ' ) ; 
snr=10/'  (snrdb/10)  ; 

%Theoretical  minimum  frequency  RMS  measurement  error 
dfrms=sqrt (3 /snr) / (pi*T) ; 
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%Generate  loop  for  sampling  trajectory  points 

%Read  positions  and  velocities  at  each  sample  point  (every  .01  sec) 
for  ts=0:T:Tmax; 
sample=ts/T+l ; 

sample=round ( sample) ;  %Sample  must  be  an  integer  value 

x=xms  (1,  sample)  ,-y=yms  (1,  sample)  ;  z  =  zms  (1,  sample)  ; 

xl=xms ( 1 , sample+1 ) ;yl=yms (1, sample+1)  ;  zl  =  zms ( 1, sample+1) ; 

v= [vxs (1, sample)  vys (1 , sample)  vzs (1, sample) ] ; 

vl= [vxs (1, sample+1)  vys (1, sample+1)  vzs (1 , sample+1) ] ; 

%Find  source  position  at  present  time 
source= [x  y  z]  ; 

%Find  source  position  at  next  time 
sourcel=[xl  yl  zl] ; 

%Find  source-to-base  station  vectors  at  present  time 
Bl=source-RL ( 1 , : ) ; 

B2=source-RL (2 , :) ; 

B3=source-RL (3 , : ) ; 

B4=source-RL (4 , :) ; 

B5=source-RL (5, :) ; 

B6=source-RL(6, : ) ; 

B7=source-RL (7 , :) ; 

B8=source-RL ( 8 , : ) ; 

B9=source-RL ( 9 ,  : )  ; 

%Find  Doppler  frequency  vector  at  present  time 
F=v/wl ; 

%Find  Doppler  frequencies  at  receivers  at  present  time  (with  noise) 
fdl=dot (B1,F) /norm(Bl) +dfrms*randn; 
fd2=dot (B2,F) /norm(B2 ) +dfrms*randn; 
fd3=dot (B3 , F) /norm(B3 ) +dfrms*randn; 
fd4=dot (B4 , F) /norm(B4) +dfrms*randn; 
fd5=dot (B5 , F) /norm(B5 ) +dfrms*randn; 
fd6=dot (B6,F) /norm(B6) +dfrms*randn; 
fd7=dot (B7 , F) /norm(B7) +dfrms*randn; 
f d8=dot (B8 , F) /norm(B8 ) +df rms*randn; 
f d9=dot (B9 , F) /norm(B9 ) +dfrms*randn; 

%Find  differential  RDOAs  at  present  time 
drdoa2a=-wl*T* (fdl-fd2)  ; 
drdoa3a=-wl*T* (fdl-fd3)  ; 
drdoa4a=-wl*T* (fdl-fd4)  ; 
drdoa5a=-wl*T* (fdl-fd5)  ; 
drdoa6a=-wl*T* (fdl-fd6)  ; 
drdoa7a=-wl*T* (fdl-fd7) ; 
drdoa8a=-wl*T* (fdl-fd8)  ; 
drdoa9a=-wl*T* (fdl-fd9) ; 

%Find  source-to-receiver  vectors  at  next  time 
Bl=sourcel-RL (1, : ) ; 

B2=sourcel-RL (2 , : ) ; 

B3=sourcel-RL (3,:); 

B4=sourcel-RL(4, : ) ; 

B5=sourcel-RL(5, : ) ; 
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B6=sourcel-RL(6,  :  )  ; 

B7=sourcel-RL(7,  : )  ; 

B8=sourcel-RL (8,  : )  ; 

B9=sourcel-RL (9, : ) ; 

%Find  Doppler  frequency  vector  at  next  time 
Fl=vl/wl ; 

%Find  Doppler  frequencies  at  receivers  at  next  time  (with  noise) 
fdl=dot (Bl, FI) /norm(Bl) +dfrms*randn; 
fd2=dot (B2,F1) /norm(B2) +dfrms*randn; 
fd3=dot (B3 , FI) /norm(B3 ) +dfrms*randn; 
f d4=dot ( B4 , FI ) /norm  ( B4 ) +df rms  *randn; 
fd5=dot (B5 , FI) /norm(B5) +dfrms*randn; 
fd6=dot (B6,F1) /norm  (B6 ) +df rms* randn; 
fd7=dot (B7 , FI) /norm(B7) +df rms* randn; 
fd8=dot (B8, FI) /norm(B8)+dfrms*randn; 
fd9=dot (B9, FI) /norm (B9 ) +df rms *randn; 

%Find  differential  RDOAs  at  next  time 
drdoa2b=-wl*T* (fdl-fd2) ; 
drdoa3b=-wl*T* (fdl-fd3) ; 
drdoa4b=-wl*T* (fdl-fd4)  ; 
drdoa5b=-wl*T* (fdl-fd5) ; 
drdoa6b=-wl*T* (fdl-fd6)  ; 
drdoa7b=-wl*T* (fdl-fd7)  ; 
drdoa8b=-wl*T* (fdl-fd8)  ; 
drdoa9b=-wl*T* (fdl-fd9) ; 

%Find  average  differential  RDOAs 
drdoa2= (drdoa2a+drdoa2b) /2  ; 
drdoa3= (drdoa3a+drdoa3b) / 2  ; 
drdoa4= (drdoa4a+drdoa4b) /2; 
drdoa5= (drdoa5a+drdoa5b) / 2; 
drdoa6= (drdoa6a+drdoa6b) /2; 
drdoa7= (drdoa7a+drdoa7b) / 2; 
drdoa8= (drdoa8a+drdoa8b) / 2; 
drdoa9= (drdoa9a+drdoa9b) / 2; 

%Find  RDOAs 

da ( 1 ) =d ( 1 ) +drdoa2 ; 
da ( 2 ) =d ( 2 ) +drdoa3 ; 
da ( 3 ) =d ( 3 ) +drdoa4 ; 
da ( 4 ) =d ( 4 ) +drdoa5 ; 
da ( 5 ) =d ( 5 ) +drdoa6 ; 
da ( 6 ) =d ( 6 ) +drdoa7 ; 
da ( 7 ) =d ( 7 ) +drdoa8 ; 
da ( 8 ) =d ( 8 ) +drdoa9 ; 

%Reset  range  differences 
d=da  ; 

%Find  RDOAs  from  geometry 

db(l)  =  norm(sourcel  -  RL(2,:))  -  norm(sourcel  -  RL(1,:)); 

db(2)  =  norm(sourcel  -  RL(3,:))  -  norm(sourcel  -  RL(1,:)); 

db(3 )  =  norm(sourcel  -  RL(4,:))  -  norm(sourcel  -  RL(1,:)); 

db(4)  =  norm(sourcel  -  RL(5,:))  -  normtsourcel  -  RL(1,:)); 
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db (5 )  =  norm(sourcel  -  RL(6,:))  -  norm(sourcel  -  RL(1,:)); 

db(6)  =  norm(sourcel  -  RL(7,:))  -  norm(sourcel  -  RL(1,:)); 

db(7)  =  norm(sourcel  -  RL(8,:))  -  norm(sourcel  -  RL(1,:)); 

db(8)  =  norm(sourcel  -  RL(9,:))  -  norm(sourcel  -  RL(1,:)); 

%Fill  matrices  for  time,  position,  and  velocity  components 
if  ts==0;t=ts+T;xx=x;yy=y;zz=z;vx=v(l) ;vy=v(2) ;vz=v(3) ; 
else  xx=  [xx  x]  ;yy=  [yy  y]  ;  zz=  [zz  z]  ;  t=  [t  ts+T]  ;  .  .  . 
vx= [vx  v(l) ] ; vy=  [vy  v(2) ] ; vz= [vz  v(3 ) ] ; 

end; 

%Use  Smith-Abel  algorithm  to  obtain  estimated  position  from  RDOAs 

%Apply  45  deg  3-dB  beam  width  of  receiver  antennas  for  signal  reception 

%Weights  of  1  define  receivers  used,  weights  of  0  define  those  not  used 

%Receivers  la, lb, lp, 10 , 11 , 12 , 14  and  15  within  view  at  time  start 
if  ts<2.12 

w=  [ 1 ; 1 ; 1 ; 1 ;  1 ; 1 ; 1 ; 0 ]  ; 

du= [ d ( 1 ) ; d ( 2 ) ;d(3) ;d(4) ;d(5)  ;d(6) ;d(7) ] ; 
dbn= [ db ( 1 ) ;db(2) ;db(3) ;db(4) ;db(5) ;db(6) ;db(7) ] ; 
source_est=smith_abel (RL,d,w,R) ; 
end; 

%Lose  Receiver  12  at  2.12  sec 
if  ts>=2 . 12  S=  ts<3 . 29 
w= [ 1 ;  1 ;  1 ; 1 ; 0 ; 1 ; 1 ; 0 ]  ; 
du= [ d ( 1 ) ;  d  ( 2 ) ;d(3) ;d(4) ;d(6) ;d(7) ] ; 
dbn= [ db ( 1 ) ;db(2) ;db(3) ;db(4) ;db(6) ;db(7) ] ; 
source_est=smith_abel (RL, d, w, R) ; 
end; 

%Lose  Receiver  14  at  3.29  sec 
if  ts>=3.29  S:  ts<4.49 
w= [ 1 ; 1 ; 1 ; 1 ; 0 ; 0 ; 1 ; 0 ]  ; 
du= [ d ( 1 ) ; d ( 2 ) ;d(3) ;d(4) ;d(7)  ]  ; 
dbn=[db(l)  ;db(2) ;db(3) ;db(4) ;db(7) ] ; 
source_est=smith_abel (RL, d, w, R) ; 
end; 

%Obtain  RDOA  for  added  receiver  from  estimated  position 
if  ts==4.49 

d(8) =norm(source_est  -  RL(9,:))  -  norm(source_est  -  RL(1,:)); 

d ( 8 ) =d ( 8 ) +drdoa9 ; 
end; 

%Lose  Receiver  15  and  gain  Receiver  16  at  4.49  sec 
if  ts>=4.49  &  ts<5.95 
w=  [ 1 ; 1 ; 1 ; 1 ; 0 ; 0 ; 0 ; 1 ]  ; 
du=[d(l) ; d ( 2 ) ; d ( 3 ) ;d(4) ;d(8) ] ; 
dbn= [ db ( 1 ) ;db(2) ;db(3) ;db(4) ;db(8) ] ; 
source_est=smith_abel (RL, d, w, R) ; 
end; 

%Lose  Receiver  11  at  5.95  sec 
if  ts>=5 . 95  Sc  ts<=Tmax 
w=  [ 1 ; 1 ; 1 ; 0 ;  0  ;  0  ;  0  ;  1 ]  ; 
du= [ d ( 1 ) ; d ( 2 )  ;  d ( 3 )  ;  d ( 8 )  ]  ; 
dbn=[db(l) ;db(2) ;db(3) ;db(8) ] ; 
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source__est=smith_abel (RL ,  d,  w,  R) ; 
end; 

%Estimated  RDOA  error  from  used  range  differences 
errt=norm(du  -  dbn) ; 

%Fill  matrix  for  estimated  RDOA  error 
if  ts==0 ; xet=errt ; 
else  xet=[xet  errt] ; 
end; 

%Estimated  position  error 

errp=norm(sourcel  -  source_est) ; 

%Fill  matrix  for  estimated  position  error 
if  ts==0 ; xep=errp; 
else  xep=[xep  errp] ; 
end 

%Return  to  trajectory  loop 
end; 

%Obtain  position  and  velocity  components  vectors 
pos=[t'  xx'  yy'  zz'] ; 
pos=pos ' ; 

v- [ t '  vx '  vy '  vz ' ] ; 
v= v ' ; 

%Output  results  to  screen 
merrt=mean(xet) ; 

fprintf { 'Mean  RDOA  Error  (ft)  =  %g  \n'  ,merrt) ; 
merrp=mean (xep) ; 

fprintf ( 'Mean  Position  Error  (ft)  =  %g  \n'  rmerrp) ; 

%Plot  graphs 

figure (1) 

plot ( t ,xet) ; grid; 

title ('RDOA  Error  vs  Time'); 

xlabel ( ' Time  -  sec ') ;ylabel (' RDOA  error  -  ft'); 

figure (2) 

plot(t,xep) ; grid; 

title (' Position  Error  vs  Time'); 

xlabel ('Time  -  sec ') ;ylabel (' Position  Error  -  ft'); 
figure (3) 

plot 3 (xx, yy, zz) ;grid; 
hold  on 

plot3  ( RL  (  :  ,  1 )  ,  RL  (  :  ,  2  )  ,  RL  (:,  3  ),'*')  ; 
hold  off 

title ( 'Missile  Trajectory  and  Remote  Receiver  Locations'); 
xlabel (' ft ' ) ;ylabel( 'ft' ) ; zlabel ( ' ft ' ) ; 
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t  ra  j  groupimc .  m 


%Position  from  FDOA  using  Smith-Abel  algorithm 

%Inputs  missile  trajectory  and  receiver  SNR  and  calculates 

%RDOA  from  FDOA  and  position  with  SNR  measurement  error 

%Random  white  Gaussian  noise  added  to  signal 

%Monte  Carlo  simulation,  15-40  dB,  600  iterations 

%Only  receivers  with  telemetry  in  3-dB  beam  width  utilized 

%trajgroupnmc .m 

%Load  White  Sands  missile  trajectory 
load  traj ; 

t= tra j ( 1 , : ) ; xm=tra j ( 2 , : ) ; ym=tra  j ( 3 , : ) ; zm= tra j (4 ,  : )  ; 

%Smooth  missile  trajectory 
xms=polyval (polyfit (t,xm, 14) , t) ; 
yms=polyval {polyfit ( t ,ym, 14) , t)  ; 
zms=polyval (polyfit (t , zm, 14) , t ) ; 

%0btain  polynomial  coefficients  for  trajectory 
xpc=polyfit (t,xm, 14)  ; 
ypc=polyfit { t ,ym, 14) ; 
zpc=polyfit (t , zm, 14) ; 

%Def ine  coefficients  for  velocity  (first  derivative  of  polynomials) 
i=l : 14 ; 

xvc= (15-i) . *xpc (i) ; 
yvc= (15-i) . *ypc (i) ; 
zvc= (15-i) . *zpc (i) ; 

% Evaluate  velocity  polynomials  to  obtain  velocity  component  vectors 
vxs=polyval (xvc, t ) ; 
vys=polyval (yvc , t ) ; 
vzs=polyval (zvc, t) ; 

% White  Sands  remote  receiver  positions  (feet) 

RL= [ -9243 . 4  -3085.7  -127;  -9193.6  -3072.3  -125;  -9278.7  -3201.8 

-92.1;  -8992.9  -2879.8  -109.3;  1757  430.4  40.4;  -6458.5  -2112.3 
-39.3;  -3704.2  -1214.7  -33.6;  607  -3517.8  208.3;  -5305.4  2653 
-153.6] ; 


%Initialize  range  difference  vectors 

d=zeros (8,1) ; 

da=zeros (8,1) ; 

di  =  zeros (8,1)  ; 

db=zeros (8,1) ; 


%Range  from  non-reference  receiver 
R=zeros (8,1) ; 


R ( 1 )  =  norm(RL(l, : ) 
R ( 2 )  =  norm(RL (1, : ) 
R(3)  =  norm(RL ( 1 , : ) 
R  (4 )  =  norm(RL (1,  : ) 
R(5)  =  norm(RL ( 1 , : ) 
R(6)  =  norm(RL (1, : ) 
R ( 7 )  =  norm(RL (1 ,  : ) 


RL  (2 , 

:)); 

RL  ( 3 , 

:)  )  ; 

RL  ( 4 , 

:)  )  ; 

RL  ( 5 , 

:)  )  ; 

RL  ( 6 , 

:))  ; 

RL  ( 7 , 

:  ) ) ; 

RL  ( 8 , 

:)  ) ; 

to  reference  receiver 
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R  ( 8 )  =  norm(RL(l,  : )  -RL(9,:)); 


%Telemetry  and  sampling  data 


f o=2 . 3e9 ;  %  Hz 

wl= (3e8/fo) / .3048;  %  feet 

T=le-2 ;  %  sec 

Tmax=6.41;  %  sec 


%Start  SNR  loop  for  receiver  SNR  values  (15-40  dB) 
for  snrdb=15:40; 

snr=10/v  (snrdb/10)  ; 

%Theoretical  minimum  frequency  RMS  measurement  error 
dfrms=sqrt (3/snr) / (pi*T) ; 

%Start  Monte  Carlo  loop  (600  iterations) 
for  count=l :  600 ; 

%Initial  range  differences,  with  Receiver  1  as  the  reference 
sourcei= [xms (1)  yms(l)  zms(l)]; 


di(l) 

=  norm(sourcei 

-  RL  (2 , 

)  ) 

-  norm(sourcei 

-  RL ( 1 , :)); 

di  (2 ) 

=  norm(sourcei 

-  RL(3, 

)  ) 

-  norm(sourcei 

-  RLU,  :)  )  ; 

di  (3 ) 

=  normfsourcei 

-  RL(4, 

)  ) 

-  normfsourcei 

-  RLU,  :))  ; 

di  (4) 

=  norm(sourcei 

-  RL ( 5 , 

)  ) 

-  norm(sourcei 

-  RLU,  :  )  )  ; 

di  (5 ) 

=  norm(sourcei 

-  RL ( 6 , 

)  ) 

-  norm(sourcei 

-  RL(1, : ) ) ; 

di  (6) 

=  normfsourcei 

-  RL (7 , 

)  ) 

-  normfsourcei 

-  RL(1, :)); 

di  (7 ) 

=  normfsourcei 

-  RL ( 8 , 

)  ) 

-  normfsourcei 

-  RL(1, : ) ) ; 

di  (8 ) 

=  norm(sourcei 

-  RL (9, 

)  ) 

-  normfsourcei 

-  RL(1, :) ) ; 

d=di ; 


%Generate  loop  for  sampling  trajectory  points 

%Read  positions  and  velocities  at  each  sample  point  (every  .01  sec) 
for  ts=0:T:Tmax; 
sample=ts/T+l ; 

sample=round ( sample ) ;  %Sample  must  be  an  integer  value 

x=xms  (1,  sample)  ;y=yms  (1,  sample)  ;  z  =  zms  (1,  sample)  ; 

xl=xms  (1,  sample+1)  ;yl=yms  (1,  sample+1)  ;  zl=zms  (1,  sample+1)  ; 

v=  [vxs  (1,  sample)  vys  (1,  sample)  vzs  (1,  sample)  ]  ; 

vl=  [vxs  ( 1 ,  sample+1 )  vys  (1,  sample+1)  vzs  ( 1 ,  sample+1)  ]  ; 

%Find  source  position  at  present  time 
source= [x  y  z] ; 

%Find  source  position  at  next  time 
sourcel=[xl  yl  zl] ; 

%Find  source-to-base  station  vectors  at  present  time 
Bl=source-RL {!,:); 

B2=source-RL (2 , : ) ; 

B3  =  source-RL (3 ,  :  )  ; 

B4=source-RL ( 4 , : )  ; 

B5=source-RL ( 5 , : ) ; 

B6  =  source-RL ( 6 , ;  )  ; 

B7=source-RL {!,:); 

B8=source-RL ( 8 ,  : )  ; 


93 


B9=source-RL (9, :) ; 


%Find  Doppler  frequency  vector  at  present  time 
F=v/wl ; 

%Find  Doppler  frequencies  at  receivers  at  present  time  (with  noise) 
fdl=dot (Bl, F) /norm(Bl) +dfrms*randn; 
fd2=dot (B2 , F) /norm(B2) +dfrms*randn; 
fd3=dot (B3 , F) /norm(B3) +dfrms*randn; 
fd4=dot (B4 , F) /norm(B4) +dfrms*randn; 
fd5=dot (B5 , F) /norm(B5)+dfrms*randn; 
fd6=dot(B6,F) /norm(B6) +dfrms*randn; 
fd7=dot (B7 , F) /norm(B7 ) +dfrms*randn; 
fd8=dot (B8,  F) /norm(B8) +dfrms*randn; 
fd9=dot (B9, F) /norm(B9) +dfrms*randn; 

%Find  differential  RDOAs  at  present  time 
drdoa2a=-wl*T* (fdl-fd2) ; 
drdoa3a=-wl*T* (fdl-fd3) ; 
drdoa4a=-wl*T* (fdl-fd4)  ; 
drdoa5a=-wl*T* (fdl-fd5)  ; 
drdoa6a=-wl *T* (fdl-fd6) ; 
drdoa7a=-wl*T* (fdl-fd7) ; 
drdoa8a=-wl*T* (fdl-fd8)  ; 
drdoa9a=-wl*T* (fdl-fd9) ; 

%Find  source-to-receiver  vectors  at  next  time 
Bl=sourcel-RL(l, : ) ; 

B2=sourcel-RL(2, : ) ; 

B3=sourcel-RL(3 , : ) ; 

B4=sourcel-RL (4 , : ) ; 

B5=sourcel-RL (5 , : ) ; 

B6=SOurcel-RL(6, : ) ; 

B7=sourcel-RL (7 , : ) ; 

B8=sourcel-RL (8 , : ) ; 

B9=sourcel-RL (9 , : ) ; 

%Find  Doppler  frequency  vector  at  next  time 
Fl=vl/wl; 

%Find  Doppler  frequencies  at  receivers  at  next  time  (with  noise) 
fdl=dot (Bl, FI) /norm(Bl) +dfrms*randn; 
fd2=dot (B2,F1) /norm(B2) +dfrms*randn; 
fd3=dot (B3 , FI) /norm(B3) +dfrms*randn; 
fd4=dot (B4 , Fl) /norm(B4) +dfrms*randn; 
fd5=dot (B5 , Fl) /norm(B5) +dfrms*randn; 
fd6=dot (B6, Fl) /norm(B6) +dfrms*randn; 
fd7=dot (B7 , Fl) /norm(B7) +dfrms*randn; 
fd8=dot (B8,F1) /norm(B8) +dfrms*randn; 
fd9=dot (B9 , Fl) /norm(B9) +dfrms*randn; 

%Find  differential  RDOAs  at  next  time 
drdoa2b=-wl*T* (fdl-fd2) ; 
drdoa3b=-wl*T* (fdl-fd3) ; 
drdoa4b=-wl*T* (fdl-fd4) ; 
drdoa5b=-wl*T* (fdl-fd5) ; 
drdoa6b=-wl*T* (fdl-fd6) ; 
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drdoa7b=-wl*T* (fdl-fd7) ; 
drdoa8b=-wl*T* (fdl-fd8)  ; 
drdoa9b=-wl*T* (fdl-fd9) ; 

%Find  average  differential  RDOAs 

drdoa2= (drdoa2a+drdoa2b) / 2; 
drdoa3= (drdoa3a+drdoa3b) / 2 ; 
drdoa4= (drdoa4a+drdoa4b) / 2 ; 
drdoa5= (drdoa5a+drdoa5b) /2; 
drdoa6= (drdoa6a+drdoa6b) / 2; 
drdoa7= (drdoa7a+drdoa7b) /2; 
drdoa8= (drdoa8a+drdoa8b) / 2; 
drdoa9= (drdoa9a+drdoa9b) /2; 


%Find  RDOAs 

da(l) 
da  (2) 
da  (3) 
da  (4) 
da  (5) 
da  { 6 ) 
da  {7 ) 
da  (8) 


=d(l) +drdoa2, 
=d(2) +drdoa3, 
=d(3) +drdoa4, 
=d(4) +drdoa5, 
=d(5) +drdoa6 
=d ( 6 ) +drdoa7 
=d{7) +drdoa8 
=d(8) +drdoa9 


%Reset  range  differences 
d=da  ; 


%Find  RDOAs  from  geometry 


db(l) 

= 

normfsourcel 

-  RL (2  , 

:)  ) 

db  (2 ) 

= 

norm(sourcel 

-  RL  (3 

:)  ) 

db  (3 ) 

= 

norm(sourcel 

-  RL(4, 

:)  ) 

db  ( 4) 

= 

normfsourcel 

-  RL ( 5 , 

:)  ) 

db  ( 5 ) 

= 

norm(sourcel 

-  RL ( 6 , 

:)  ) 

db  ( 6 ) 

= 

norm{sourcel 

-  RL(7, 

:)  ) 

db  ( 7 ) 

= 

normfsourcel 

-  RL ( 8 , 

:)  ) 

db  ( 8 ) 

= 

no:rm  ( sourcel 

-  RL ( 9 , 

:)  ) 

norm(sourcel 

norm(sourcel 

norm(sourcel 

norm(sourcel 

norm(sourcel 

norm(sourcel 

norm(sourcel 

norm(sourcel 


-  RL ( 1 ,  : )  )  ; 

-  RL ( 1 , : ) )  ; 

-  RLU,  : )  )  ; 

-  RL ( 1 ,  :  )  )  ; 

-  RL ( 1 ,  : ) )  ; 

-  RL ( 1 ,  :  )  )  ; 

-  RL ( 1 , : ) )  ; 

-  RLU,:)); 


%Use  Smith-Abel  algorithm  to  obtain  estimated  position  from  RDOAs 
% Apply  45  deg  3-dB  beam  width  of  receiver  antennas  for  signal  reception 
%Weights  of  1  define  receivers  used,  weights  of  0  define  those  not  used 
%Receivers  la, lb, IP/ 10 , 11, 12 , 14  and  15  within  view  at  time  start 
if  ts<2 . 12 

w= [ 1 ; 1 ; 1 ; 1 ; 1 ; 1 ; 1 ; 0 ] ; 

du= [ d ( 1 ) ; d ( 2 ) ; d ( 3 ) ; d ( 4 ) ; d ( 5 ) ; d ( 6 ) ; d ( 7 ) ] ; 
dbn=  [  db  ( 1 )  ;  db  ( 2 )  ;  db  ( 3 )  ;  db  ( 4 )  ;  db  ( 5 )  ;  db  ( 6 )  ;  db  ( 7  )  ]  ; 
source_est=smith_abel  (RL,  d,  w,  R)  ; 
end; 


%Lose  Receiver  12  at  2 . 12  sec 
if  ts>=2 . 12  &  ts<3 . 29 
w= [ 1 ; 1 ; 1 ; 1 ; 0 ; 1 ; 1 ; 0 ] ; 
du= [ d ( 1 ) ; d ( 2 ) ; d ( 3 ) ; d ( 4 ) ; d ( 6 ) ; d ( 7 )  ]  ; 
dbn=  [ db  ( 1 )  ;db(2);db(3);db(4);db(6);db(7)]; 
source_est=smith_abel  (RL,d,w,R)  ; 
end ; 
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%Lose  Receiver  14  at  3.29  sec 
if  ts>=3.29  Sc  ts<4.49 
w=  [ 1 ; 1 ; 1 ; 1 ; 0 ;  0 ;  1 ;  0  ]  ; 
du= [ d { 1 ) ; d ( 2 ) ; d ( 3 ) ; d { 4 ) ; d ( 7 ) ] ; 
dbn=[db(l) ;db(2);db(3);db(4);db(7)]; 
source_.es t =smi th_abel  (RL,  d,w,R) ; 
end; 

%Obtain  RDOA  for  added  receiver  from  estimated  position 
if  ts==4.49 

d  ( 8 )  =norm  ( source_est-RL  ( 9 ,  : )  )  -norm  { source_est-RL  ( 1 ,  *. )  ) 
d ( 8 ) =d ( 8 ) +drdoa9 ; 
end; 

%Lose  Receiver  15  and  gain  Receiver  16  at  4.49  sec 
if  ts>=4.49  Sc  ts<5.95 
w= [ 1 ; 1 ; 1 ; 1 ; 0 ; 0 ; 0 ; 1 ] ; 
du= [ d ( 1 ) ; d { 2 ) ; d ( 3 ) ; d ( 4 ) ; d ( 8 ) ] ; 
dbn= [ db ( 1 } ;db{2) ;db{3) ;db(4) ;db(8) ] ; 
source_est=smith_abel  (RL,  d,  w,  R)  ; 
end; 

%Lose  Receiver  11  at  5.95  sec 

if  ts>=5.95  Sc  ts<=Tmax 
w= [ 1 ; 1 ; 1 ; 0 ; 0 ; 0 ; 0 ;  1  ]  ; 
du=[d(l) ; d ( 2 ) ; d { 3 ) ;d(8) ]  ; 
dbn= [db ( 1 ) ;db(2);db(3);db(8)]; 
source_est=smith_abel (RL, d# w, R) ; 
end; 

%Estimated  RDOA  error  from  used  range  differences 
errt=norm(du  -  dbn) ; 

%Fill  matrix  for  estimated  RDOA  error 
if  ts==0;xet=errt; 
else  xet=[xet  errt] ; 
end; 

%Estimated  position  error 

errp=norm (source 1  -  source_est) ; 

%Fill  matrix  for  estimated  position  error 
if  ts==0;xep=errp; 
else  xep=[xep  errp] ; 
end 

%Return  to  trajectory  loop 
end; 

%Calculate  mean  RDOA  and  position  errors  for  missile  trajectory 

%Fill  matrices  for  these  values 
merrt=mean (xet ) ; 
merrp=mean (xep) ; 

if  count==l;meanxet=merrt ;meanxep=merrp; 

else  meanxet=  [meanxet  merrt]  ;meanxep=  [meanxep  merrp];end 
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%Return  to  Monte  Carlo  loop 
end; 

%Calculate  and  print  mean  values  for  Monte  Carlo  iterations 
mcet=mean (meanxet ) ; 
mce=mean(meanxep)  ; 
f print f ( ' SNR(dB)  =  %g\n' ,  snrdb); 
fprintf (' Error  in  RDOA  (ft)  =  %g\n',  meet); 
fprintf ( 'Error  in  Position  (ft)  =  %g\n' ,  mce) ; 

%Fill  SNR,  RDOA,  and  position  error  matrices  for  plotting 
if  snrdb==15 ; snrdbm=snrdb;mcetm==mcet ;mcem=mce; 

else  snrdbm=  [  snrdbm  snrdb]  ;mcetm=  [mcetm  meet]  ;mcem=  [mcem  mce];  end; 

%Return  to  SNR  loop 
end; 

%Save  data  to  workspace 

save  snrdata2  snrdbm  mcetm  mcem; 

%Plot  graphs 
figure (1) 

plot (snrdbm, mcetm) ;grid; 

title ('Range  Difference  Error  vs.  SNR'); 

xlabel ( ' SNR  -  dB' ) ;ylabel ( 'Range  Difference  Error  -  ft'); 
figure ( 2 ) 

plot (snrdbm, mcem) ;grid; 

title ( 'Position  Error  vs.  SNR'); 

xlabel ('SNR  -  dB ') ;ylabel (' Position  Error  -  ft'); 
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APPENDIX  B.  RECEIVER  AZIMUTH  AND  ELEVATION  ANGLES 


Figures  B.l  -  B.9  show  the  azimuth  (theta)  and  elevation  (phi)  angles  vs.  time  for 
each  of  the  nine  receivers  on  the  White  Sands  Missile  Range.  These  plots  were  obtained 
from  the  MATLAB  program  angles.m  whose  source  code  is  listed  in  Appendix  A. 

Angles  for  Receiver  la  -  Reference  Receiver 


Figure  B.l.  Reference  Receiver  Azimuth  and  Elevation  Angles 
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Angles  for  Receiusr  1 4 


Figure  B.6.  Receiver  12  Azimuth  and  Elevation  Angles 


Angles  for  Receiver  1 2 


Figure  B.7.  Receiver  14  Azimuth  and  Elevation  Angles 


Time  -  sec 


Figure  B.8.  Receiver  15  Azimuth  and  Elevation  Angles 

Angles  for  Receiver  16 


1  2  3  4  5  6 

Time  -  sec 

Figure  B.9.  Receiver  16  Azimuth  and  Elevation  Angles 
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